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We show that dipolar interactions between ultracold polar alkali dimers in optical lattices can 
be used to realize a highly tunable generalization of the t-J model, which we refer to as the t-J- 
V-W model. The model features long-range spin-spin interactions J z and J± of XXZ type, long- 
range density-density interaction V, and long-range density-spin interaction W, all of which can be 
controlled in both magnitude and sign independently of each other and of the tunneling t. The "spin" 
is encoded in the rotational degree of freedom of the molecules, while the interactions are controlled 
by applied static electric and continuous-wave microwave fields. Furthermore, we show that nuclear 
spins of the molecules can be used to implement an additional (orbital) degree of freedom that is 
coupled to the original rotational degree of freedom in a tunable way. The presented system is 
expected to exhibit exotic physics and to provide insights into strongly correlated phenomena in 
condensed matter systems. Realistic experimental imperfections are discussed. 

PACS numbers: 67.85.-d, 71.10.Fd, 33.80.-b, 33.20.-t 



I. INTRODUCTION 

Ultracold diatomic polar molecules have recently at- 
tracted a great deal of attention both experimentally 
and theoretically [1-6]. Two features of diatomic po- 
lar molecules make them particularly interesting as com- 
pared to the more typical systems of ultracold alkali 
atoms. First, polar molecules possess a permanent dipole 
moment, which can be manipulated with external fields 
and which can lead to long-range anisotropic interac- 
tions. This contrasts with atoms whose interactions are 
typically short-range and isotropic. Second, the inter- 
nal level structure of diatomic polar molecules is much 
richer than that of atoms and, although more difficult 
to control, allows, in principle, for richer physics. These 
two features make diatomic polar molecules attractive 
for numerous applications including quantum computa- 
tion, quantum simulation, precision measurements, and 
controlled quantum chemistry [1-5]. 

In Ref. [7], it was shown that these two unique fea- 
tures allow ultracold polar molecules in optical lattices 
to simulate a highly tunable generalization of the t-J 
model [8-10] referred to as the t-J-V-W model. This 
model makes use of the rotational degree of freedom of 
the molecules and features tunneling t, density-density 
interaction V, density-spin interaction W , and spin-spin 
interactions J z and J± of XXZ type. As a first step to- 
wards understanding this model, Ref. [7] showed that the 
simplest experimentally realizable case of the t-J-V-W 
model with V — W — J z — allows to strongly enhance 
the superconducting (i.e. superfluid for our neutral sys- 
tem) region of the ID phase diagram relative to the usual 
t-J model. 

In the present paper, we provide the details behind 
the derivation of the t-J-V-W model. In particular, we 
show that the manipulation of the rotational degree of 
freedom of the molecules via DC electric and microwave 



fields allows to achieve full control of the coefficients of 
the t-J-V-W Hamiltonian and discuss the implications of 
this control on the accessible manybody physics. Specifi- 
cally, one can tune the system into exhibiting the physics 
very similar to the original t-J model, whose phase dia- 
gram is still highly controversial beyond one dimension 
[8-10]. Alternatively, one can access a wide range of other 
regimes that include the spin-1/2 XXZ magnet and nu- 
merous extensions of the t-J model, some of which are 
believed to exhibit enhanced superfluid correlations. We 
also show how to control the spatial anisotropy of the 
Hamiltonian by changing the direction of the applied DC 
electric field and how to control the optical potential ex- 
perienced by different rotational states by an appropriate 
choice of lattice laser beams. 

Furthermore, we study in detail the generalization of 
the t-J-V-W model to the case where not only the rota- 
tional degree of freedom of the molecules, but also their 
nuclear degrees of freedom play an important role. Due 
to the relative simplicity of their production, the only 
ultracod polar molecules currently available in their elec- 
tronic, vibrational, and rotational ground states are al- 
kali dimers KRb [11-14] and LiCs [15, 16]. Therefore, we 
focus on the hyperfine structure of alkali dimers, which 
has recently been studied theoretically [17-21] and ex- 
perimentally [22-24]. Specifically, we show how the ap- 
plied DC electric field can be used to couple and decouple 
rotational and nuclear degrees of freedom, thus allow- 
ing for the control of nuclear spin effects. In the case 
where nuclear spins are coupled to the rotational degree 
of freedom, we show that the nuclear spins can function 
either as classical - possibly spatially dependent - mag- 
netic field or as a separate (orbital) quantum degree of 
freedom with a highly tunable interaction with the rotor. 
We also point out possible promising applications of the 
system to quantum information processing. Since ultra- 
cold ground-state polar alkali dimers are already avail- 
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able in experiments [11-16] and are even loaded in opti- 
cal lattices [13], we expect our results to be immediately 
applicable to current experiments. 

Our work builds on an extensive body of literature 
studying the many-body dynamics of polar molecules in 
a lattice and making use of the internal rotational struc- 
ture [3, 21, 25-36] and of fine and hyperfine structure of 
molecules with a single electron outside a closed shell 
[26, 27, 33]. We would also like to specifically high- 
light recent work in Refs. [21, 36], which make important 
steps towards the understanding of the effects of hyper- 
fine structure on many-body physics with alkali dimers. 

The remainder of the paper is organized as follows. In 
Sec. II, we introduce the t-J-V-W Hamiltonian in the 
presence of both a rotational and a nuclear degree of 
freedom and describe its main features. Then, in Sees. 
Ill - VI, we present a detailed derivation and discussion 
of this Hamiltonian. In particular, in Sec. Ill, we study 
the rotational and hyperfine structure of the molecules in 
the presence of a DC electric field. In Sec. IV, we study 
the optical potential and the associated tensor shifts. In 
Sec. V, we use the results of Sees. Ill and IV to give 
a detailed derivation of the final Hamiltonian. In Sec. 
VI, we find the regimes, in which the model is stable 
to loss via chemical reactions. Finally, in Sec. VII, we 
present the conclusions. Appendix A presents formulas 
useful for studying the single-molecule Hamiltonian and 
dipole-dipole interactions between molecules. Appendix 
B describes the phenomenon of interaction-assisted tun- 
neling, which arises if one considers small corrections to 
the t-J-V-W model. 



II. THE HAMILTONIAN AND ITS FEATURES 

In this Section, we introduce the t-J-V-W Hamiltonian 
in the presence of both a rotational and a nuclear degree 
of freedom and describe its main features. The detailed 
derivation is postponed until Sees. Ill- VI. 

We consider diatomic polar molecules confined to a 
single plane (e.g. using a strong ID optical lattice) and 
subject to a DC electric and, possibly, one or more 
continuous-wave (CW) microwave fields. Furthermore, 
in that plane, the molecules are assumed to be loaded in 
the lowest band of a 2D optical lattice. Such a system 
is not far out of reach experimentally: indeed, loading 
of KRb molecules into ID [13] and 3D [37] lattices and 
of homonuclear Cs2 molecules into a 3D lattice [24] has 
already been demonstrated. 

As shown in Sees. III-V, taking into account the ap- 
plied DC and microwave fields, we can reduce the internal 
structure of each molecule to a tensor product of a two- 
level dressed rotational degree of freedom (dressed states 
labeled by |mo) and |mi); angular momentum operator 
on site j labeled by Sj ) and a two- level nuclear degree of 
freedom (states labeled by | f) and | J,) ; angular momen- 
tum operator on site j labeled by Tj). In Sees. Ill- VI, 



we derive the following Hamiltonian: 



(i,j)ma 



-j5> dd (R 2 -R,) 



+Vn inj +W{n l S z J + njS*) 



h.c. 



J z Sf t S Z j + — (S^Sj + Si Sj~) 



+ay,s*t:. (i) 



This Hamiltonian, together with the full control over its 
coefficients and with the detailed study of the hyperfine 
structure, is the main result of the present paper. The 
first term (oc t m ) describes tunneling of molecules; the 
second term (oc V^d) describes dipole-dipole interactions; 
while the last term (oc A) describes hyperfine interac- 
tions. Let us describe each of these terms, including the 
necessary definitions and the physical origin. 

Let us begin with the tunneling term oc t m . The 
bosonic or fermionic creation operator cj- mcr creates a 
molecule on site j in the dressed rotor state to (= too 
or mi) and nuclear state a (=f or |). The notation 
indicates the sum over nearest neighbors, where each pair 
of nearest neighbors is included only once. Throughout 
the paper, we set H = 1. As we will show in Sec. IV, 
the tunneling amplitudes t mo and t mi can be made ei- 
ther equal or different by choosing the polarization and 
frequency of the optical fields creating the lattice and by 
choosing the dressed states |too) and |mi). The over- 
all magnitude of the amplitudes t m can be tuned from 
zero up to a few kHz by changing the intensity of the 
optical fields. Notice that the chemical potential is not 
included in Eq. (1) since H conserves the total number 
of molecules in each internal state |toct) and since the 
system is not coupled to a reservoir of molecules. The 
effect of a chemical potential can be modeled by control- 
ling the total number of molecules in each internal state 
\ma) during the preparation stage. 

Let us now describe the dipole-dipole interaction term 
oc Vdd- The operator nj ma = Cj ma Cj mrT counts the num- 
ber of molecules on site j in the dressed rotor state m 
and nuclear state a, while the operator n jm = rij m<7 
counts the number of molecules on site j in the dressed 
rotor state to irrespective of the nuclear state. We assume 
[see Sec. VI] that on-site interactions and/or on-site de- 
cay for two molecules are so large that molecules obey 
the hardcore constraint, i.e. each site can be occupied 
by either or 1 molecules [although it is straightforward 
to extend the model to finite on-site interactions (see 
e.g. Ref. [25])]. The operators — (rij mn — 



n jmi )/2, 

S t = E/jmo/imi^ and Sj = (5t)t are the usual 
spin- 1/2 angular momentum operators on site j describ- 
ing the two-level dressed rotor degree freedom and satis- 
fymg[SJSf]=±Sf. 

As shown in Fig. 1, the 2D plane, which the molecules 
are confined to, is assumed to be the X-Y plane, while 
the vector perpendicular to it defines the Z-axis (note 
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FIG. 1. (color online). The geometry of the setup. The 
molecules are assumed to be in the X-Y plane. A typical 
vector R in that plane has polar coordinates (R, <E>). The 
direction of the DC electric field has spherical coordinates 
(0q, 5>q) in the X-Y-Z coordinate system. The quantization 
axis z for the spins lies along the applied DC electric field. 
The other two axes (x and y) of the spin coordinate system 
are not shown. The cosine of the angle between R and i is 
equal to R ■ z = sinOo cos($ — $o) [see Eq. (2)]. 



the use of the upper case to denote the spatial axes). 
All angular momenta are, on the other hand, quantized 
along the z-axis (note the use of the lower case to denote 
angular momentum axes), which is the axis along which 
the DC electric field is applied. The z-axis has spherical 
coordinates (Oo^o) relative to the X-Y-Z coordinate 
system. In Eq. (1), R^ is the position of site i in the X-Y 
plane. Classical dipole-dipole interaction energy between 
two unit electric dipoles oriented along z and located at 
sites i and j is then given by 



^id(R) 
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[l-3sin 2 e cos 2 ($-$ )] , (2) 
(R, <f>) in polar coordinates, and 



where R = Rj — R 3 
R = R/i? is a unit vector along R. In Eq. (2), R • z = 
sin Go cos($ — $o) is the cosine of the angle between R 
and z. Vdd(R) is used in Eq. (1) and reproduces the 
usual dipole-dipole interaction behavior with head-to-tail 
attraction when z = R and side-to-side repulsion when 
z = Z. The dipole-dipole interaction term in Eq. (1) is 
multiplied by 1/2 since we double-count. 

The origin of the J z , J±, V, and W terms in Eq. (1) 
can be evinced with the following simple example, which 
does not involve the application of microwave fields. The 
rotational degree of freedom of a single molecule is de- 
scribed by the angular momentum operator N. Let us 
pick as | mo) and \m-\) the lowest two N z — states of 
the molecule in the presence of a DC electric field along 
z. Due to the applied electric field, these states are not 
eigenstates of N 2 and possess nonzero dipole moments. 
One can then intuitively think of the ground state |mo) 
as a dipole /xq = hqz oriented along the DC field (i.e. 



> 0) and of the excited state | mi) as a dipole n\ = fiiz 
oriented against the DC electric field (i.e. \i\ < 0). Let 
us now consider classical dipole-dipole interaction energy 
Edd between a dipole fii — (/io n mi + l jL i n im 1 ) z a t site i 
and a dipole fij = (/z*mjm + Minimi) 2 at site j, where 
nkm indicates whether the molecule on site k is in state 
\m) (n km = 1) or not (n km = 0): 

Edd = 4-TreolR 1 - R-| 3 ^ ' ^ ~ ^ ' ^ H ' ^ ® 

= Vrfd(Ri - Rj)(^o«im + Minimi )(A t 0"imo + 

= V dd (Ri - Rj) [J Z S*S* + Vmrij + W( ni S* + njSf)] , 

where J z = (/z - Mi) 2 , V = (fMi + Mi) 2 /4, and W = 
(/io — Mi)/2. The V term describes density-density inter- 
actions, and is the only term that survives if one averages 
Edd over the internal states |mo) and |mi) of each of the 
two molecules. Furthermore, only the V term survives 
if /iq — Hi, in which case dipole-dipole interaction can- 
not depend on the internal states of the two molecules. 
The J z term describes an Ising-type spin-spin interac- 
tion. Since J z is non-negative in this example, it favors, 
for V dd > 0, antialignment of molecules on sites i and 
j. This makes sense since two side- by-side dipoles repel 
if they are aligned, but attract if they are antialigned. 
Finally, the W term describes spin-density interaction. 
In the language of quantum magnetism, the presence of 
a molecule on site i creates, via the term WriiS^, an 
effective magnetic field along z for the spin on site j. 
As one can see from this discussion, an important differ- 
ence of our Hamiltonian from Refs. [26, 27], which also 
engineer magnetic models using molecules in optical lat- 
tices, is that we use dipole-dipole interactions in first or- 
der (rather than second order), which allows for stronger 
interactions. 

To understand the origin of the Jj_ term in Eq. (1), 
one has to take into account the transition dipole mo- 
ment /ioi between |mo) and |mi) [25]. In the same 
way, in which an optical excitation can be exchanged 
between two two-level atoms that are within an opti- 
cal wavelength of each other [38], the J± term describes 
the exchange of a microwave excitation; in this exam- 
ple, J± = 2/Iqi- For a molecule on site i and a molecule 
on site j that share one microwave excitation, the Jj_ 
term is diagonalized by the symmetric and antisymmet- 
ric states (|momi)jj ±1 mimo)y)/ v2, which would be the 
microwave equivalent of optical superradiant and subra- 
diant states [38]. The presence of the J± term is one of 
the main differences of Eq. (1) from the Hamiltonian dis- 
cussed in Ref. [21]. While this simple example illustrates 
the physical origin of the dipole-dipole interaction terms 
featured in Eq. (1), we will show in Sec. VB, that this 
form of dipole-dipole interactions is much more general. 
In particular, we will show that it applies even when mi- 
crowave fields are applied and when states with N z ^ 
are involved. 

The recipe for controlling dipole-dipole interactions 
with applied DC electric and microwave fields is one of 
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the main results of the present paper. Specifically, as we 
will discuss in Sec. II A 1, the spatial anisotropy of the 
interactions can be controlled via (0 O , ^o)- More impor- 
tantly, as we will discuss in Sees. II A 1, II A 2 and derive 
in Sec. V B, by tuning the strength of the DC electric field 
as well as the frequency and intensity of the applied mi- 
crowave field(s), one can achieve complete control over 
signs and relative amplitudes of the coefficients V, W, 
J z , and Jj_. The strength of the resulting dipole-dipole 
interactions (quoted for nearest neighbors separated by 
500 nm) is ~ 0.4 kHz in KRb and - 40 kHz in LiCs. 
These dipole-dipole interactions (particularly in the case 
of LiCs) are substantially stronger than superexchangc 
interactions in cold atoms (<S! 1 kHz [39]), making mag- 
netism easier to access in such molecular systems than in 
atomic systems. 

Finally, let us describe the hyperfine interaction term 
(oc A) in Eq. (1). The operator n J(T = £ rij ma counts 
the number of molecules on site j with nuclear spin a 
irrespective of the rotational state. The operators T 2 = 



(n it - nji )/2, T+ = £ 



m jm 



and T; 



are the usual spin- 1/2 angular momentum operators on 
site j describing the two-level nuclear degree freedom and 
satisfying \r?,T~p\ — ±Tp. The hyperfine interaction of 
the form ASfT? relies on the fact [see Sec. Ill] that, for 
a generic DC electric field, the hyperfine interaction can 
be projected on states |mo) and \m\) and, moreover, is 
diagonal in the same nuclear spin basis in both states (the 
basis, in which the two nuclei are decoupled from each 
other). Thus ASfT? simply reflects the fact that the 
energy difference between any two of these eigenstates 
(| "f) and | I)) is generally not the same in \mo) and 
| mi): the flip of the nuclear degree of freedom from | J,) 
do | f) in l m o) takes an energy larger by an amount A 
than in \mi). In Sec. VA, we show that the hyperfine 
interaction constant A can be tuned, via the strength of 
the DC electric field and via the choice of nuclear spin 
states, from zero to almost any value up to ~ 1 MHz in 
KRb and up to ~ 100 kHz in LiCs. Moreover, as we will 
note in Sees. II B and V A, while the interaction S*T? is 
the easiest form of the hyperfine interaction that one can 
obtain, any interaction between Si and Tj is, in principle, 
achievable. 

Both fermionic ( 40 K 87 Rb [11]) and bosonic ( 7 Li 133 Cs 
[15] and K Rb [14]) species are available experimen- 
tally. The bosonic [21, 25, 30-32, 40-49] and fermionic 
[21, 29-32, 36, 48, 50, 51] cases are expected to give rise 
to different physics. 

If tunneling in the third direction is negligible or if 
stabilization against collapse and/or chemical reactions 
in the third direction can be achieved without strong 
dipole-dipole repulsion (see Sec. VI), we can extend the 
Hamiltonian to 3D. A ID geometry [7], as well as non- 
square lattices can also be considered. 



A. Rotational degree of freedom alone 

In this Section, we ignore the nuclear degree of free- 
dom in Eq. (1) and discuss the tunability of the resulting 
model, as well as the physics that can be accessed with 
it. 



1. Quantum magnetism 

In this Section, we further suppose that the tunneling 
is negligible. The simplest scenario is then the case of a 
single molecule per site. In this case, rij = 1 for all sites 
i. This means that the term in the Hamiltonian pro- 
portional to V is a constant and can be dropped. The 
term proportional to W gives an effective magnetic field 
on each site. Ignoring edge effects, this magnetic field is 
uniform, making the W term commute with the Hamil- 
tonian. In this case, the W term can also be ignored, so 
that Eq. (1) reduces to 



H - - ^ v dd (R< - Rj 



J 



■ L fQ+ 



-S-Sf)U) 



The important features of the interaction in Eq. (4) are 
that it is long-range, anisotropic in both space and spin, 
and highly tunable via the magnitude of the DC elec- 
tric field, O , 7 the choice of rotational states, and the 
number, frequency, and intensity of applied microwave 
fields. In Ref. [52] and Refs. [31, 32], this Hamiltonian 
is studied in the ID geometry in the context of ions and 
molecules, respectively. The J z = case is also studied 
in the context of molecules in Ref. [25] . In Ref. [35] , this 
Hamiltonian is studied in the context of exciton-impurity 
interactions generated with polar molecules. Related lat- 
tice models with dipolar interactions are also studied in 
the context of Frenkel excitatons [53, 54]. In Ref. [55], a 
similar Hamiltonian is studied in the context of molecular 
Wigner crystals for quantum memory applications. 

As we will show in Sec. VB, if we parametrize J z and 
J_L as J z — |J|cos-0 and J± = |J|sim/>, any value of 
tp can be achieved by an appropriate combination of DC 
electric and microwave fields. In other words, one can ac- 
cess the full parameter space. For example, one can get a 
classical Ising model with J± — 0, a pure XX model with 
J z = 0, or the SU(2)-symmetric Heisenberg interaction 

(Jz = J±). 

By changing the direction (Go, <I>o) of the applied elec- 
tric field, one can control the spatial anisotropy of the 
interaction [25, 56, 57]. In particular, one can set to 
zero couplings along one or two directions in the X-Y 
plane. We assume a 2D square-lattice geometry and de- 
fine V x ,y = [1-3 sin 2 e cos 2 ($ Xi y-$ )](A 2 + y 2 )- 3 / 2 
as the coupling coefficient between the origin and the 
site with coordinates (X, Y) . The coordinates are given 
in units of lattice spacing a = A/2, where A is the 
wavelength of the light used to form the lattice. Here 
<&x.y = Arg(A + iY) is the polar angle of the vector 



(X,Y) in the plane, and (0o,"I>o) ar e the polar and az- 
imuthal angles of the applied DC electric field in the 
(X, Y, Z) coordinate system. 

In three dimensions, two cones making an angle 
cos _1 (l/\/3) ~ 0.30-7T with the applied DC electric field 
(which points along z) give vanishing dipole-dipolc inter- 
actions. As we tilt the electric field from Z towards the 
X-Y plane (i.e. increase ©o), interactions in the plane 
start changing magnitude in an anisotropic fashion. In 
particular, when sinGo = 1/V3 (Oo ~ 0.207r), the cone 
of vanishing interaction touches the X-Y plane giving one 
line in the plane along which dipole-dipolc interactions 
vanish. As shown in Figs. 2(a) and (b), using $ = or 
$o = 7r /4, we can set Vi t o = or Vx,i = 0, respectively. 
An interesting feature of setting Vi t x = is that (pro- 
vided interactions beyond Vi,_i are ignored), this turns 
a square lattice into an effective triangular lattice. We 
note that the interactions Vx,y in Fig. 2 are normalized 
by the magnitude of the largest one. 

As we tilt the electric field further, the single line 
of vanishing interaction splits into two, and the angle 
between the two lines increases up to a maximum of 
2cos~ 1 (l/v / 3) ~ 0.61-7T when z (i.e. the electric field) 
is in the X-Y plane. This way, for example, coupling 
along two orthogonal directions can be set to zero when 
sin Go = \/2/3. In particular, as shown in Fig. 2(c), at 
<I>o = 0, we have Vi,i = V\-\ = 0, while Vq,i i g un- 
changed and while the sign of Vi t o flips. Alternatively, as 
shown in Fig. 2(d), at $o = 7r /4, we get Vi,o = Vb,i = 0, 
while Vi t —i is unchanged and while the sign of Vx,i is 
flipped. 

Finally, some couplings can be set equal to each other. 
For example [Fig. 2(e)], at sin0 o « 0.51 and $ = 0, 
we have V h0 = Vi,i « 0.2Vb,i. Alternatively [Fig. 2(f)], 
at sin 0o = 0.66 and $o = 7r /4, we have 0.29Vb.i = 
0.29V_i,i = ~Vi,i. 

In experiments, it is difficult to achieve a perfect occu- 
pation of exactly one molecule per site. There will, thus, 
always be empty sites, which will play the role of defects 
in the corresponding spin model. Furthermore, empty 
sites can be introduced on purpose to emulate the effect 
of static non-magnetic impurities in quantum magnets. 



2. The t-J-V-W model 

Allowing for tunneling between sites, we arrive at the 
following Hamiltonian, which we refer to as the t-J-V-W 
model: 
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FIG. 2. (color online). Control over the spatial anisotropy 
of the interactions. The numbers near particular sites show 
the dipole-dipole interaction coefficient Vx,y (normalized by 
the magnitude of the largest Vx,y) between the origin and 
the site (X,Y). (a) 9 = sin -1 (l/V3) « 0.20tt, $ = 0. 
(b) Bo = sin-^l/v^), $o = it/4, (c) 6 = ahr 1 (y/k/S) « 
0.30?r, $o = 0. (d) 9o = sin-^v^), $ = tt/4. (e) 8 = 

sm _1 (Y (2y/2 - l)/(6v/2 - 3/2)) w 0.17tt, $ = 0. (f) 6 = 
sin _1 (v/(2- l/V2)/3)) w 0.23tt, $ = tt/4. 



H = - *« 

{i,j)m<r 



&4 v-i-i fT C 7 



h.c. 



+Vn l n 1 + W(mS? + UjSf) 



(5) 



This model is an extension of the t-J model [58-60]. 
The t-J model emerges from the large-JJ expansion of 
the Hubbard model. Despite significant efforts to iden- 
tify the phase diagram of the t-J model, only the ID 
phase diagram is relatively well-established (via numer- 
ical methods) [58-60]. It has also been demonstrated 
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that, in ID, the addition of repulsive nearest-neighbor 
interactions V ^ ■ n,-nj+i (giving rise to the so-called t- 
J-V model) and next-nearest-neighbor repulsive interac- 
tions V J2i n i n i+2 (giving rise to the so-called t-J-V-V 
model) can strengthen superconducting (i.e. superfluid 
for our neutral system) correlations in the t-J model [58]. 
Ref. [58] also argues that this effect will manifest itself in 
a 2D geometry as well. A confirmation of this statement 
can have important implications in the understanding of 
high-temperature superconductivity in cuprates. 

The highly tunable model in Eq. (5) provides unique 
opportunites to study a generalized t-J model in ID and 
2D geometries. Some of the important features of Eq. 
(5) as compared to the t-J model are as follows. First, 
instead of antiferromagnetic nearest-neighbor Heisenberg 
interactions (J± = J z > 0), Eq. (5) features long-range 
(1/i? 3 ) XXZ interactions with values of J± and J z that 
can be independently tuned in magnitude and sign. For 
example, by adjusting the sign of J±, one can obtain 
the unusual ferromagnetic interactions for fermions and 
antiferromagnetic interactions for bosons. Second, in- 
stead of the (— j + V) Ylu,j) n i n j interactions, Eq. (5) 
features long-range (1/i? 3 ) density-density interactions 
(oc V), which can be easily made repulsive to favor su- 
perfluid correlations. Third, t in Eq. (5) can be tuned 
independently from J z , J±, and V, and W. In partic- 
ular, one can access the regime | J z \, \ J±\ > t, which is 
not possible if J ~ t 2 /U. Finally, the term oc W, which 
describes density-spin interactions, can be made nonzero 
and can compete with spin-spin interactions (oc Jj_, J z ) 
and thus favor new types of spin ordering. 

Parameter regimes close to the original t-J and t-J- 
V-V models can be achieved. In particular, to obtain 
the model most similar to the t-J model, we show in 
Sec. V B how to set W = and J z = J± = -W > 0. 
We also show how to set W — 0, J z = J± > 0, and 

V = 0.1J Z , which is expected to result in a suppression 
of phase separation relative to the t-J model [58]. 

Being a generalization of the already highly nontrivial 
t-J model (particularly beyond ID), the Hamiltonian in 
Eq. (5) is expected to give rise to very rich many-body 
physics. Only a narrow range of this physics has been 
studied so far. In particular, the Hamiltonians consid- 
ered in Refs. [4, 5, 40-47, 49-51] are reminiscent of the 
restriction of Eq. (5) to a single rotational state. The use 
of more than one rotational state for manybody physics 
with diatomic polar molecules has been considered be- 
fore in Refs. [3, 21, 25-36, 48, 61-64]. Finally, in Ref. 
[7], using density matrix renormalization group (DMRG) 
[65-67], we studied the ID phase diagram of the sim- 
plest experimentally realizable regime of Eq. (5), where 

V = W = J z = 0, t mo = t mi = t, and the two re- 
maining parameters are molecule density and J±/t . As 
expected from the above discussion, we indeed found an 
enhancement of superfluid correlations and a suppression 
of phase separation relative to the usual t-J model. 

Preparation of the phases can be done, for example, 
by applying an additional microwave field coupling the 



two dressed rotor states and performing an adiabatic pas- 
sage from a state that is easy to prepare to the desired 
ground state by tuning the Rabi frequency and the de- 
tuning of the microwave field [32] . This extra microwave 
field, which gives rise to terms proportional to Y] j Sj and 
J2j Sj (i- e - effective x and z magnetic fields), can also be 
thought of as a way of enabling the simulation of a richer 
class of models where J^j Sj i s not conserved. We ex- 
pect that, by analogy with Ref. [32], preparation of the 
phases of interest can often be done without single-site 
addressability. 

Molecules in the rovibrational ground state can be 
detected by converting them back to atoms [11]. Fur- 
thermore, efforts towards achieving optical cycling in 
molecules are under way [68, 69]. There is, thus, hope 
that powerful tools for the detection of molecular phases 
can be borrowed [32] from experiments with ultracold 
atoms. These tools include noise-correlations in the time- 
of-fiight absorption imaging [70-72] and direct in-situ flu- 
orescent imaging [73, 74]. In Ref. [7], the possibility of 
probing the phase diagram of Eq. (5) with center-of-mass 
Bloch oscillations [75, 76] is also discussed. 

The model can be extended to more than two dressed 
rotor states. By applying a sufficient number of mi- 
crowave fields, one can achieve significant tunability of 
the coefficients even in the resulting more complicated 
models. 



B. Effects of the nuclear degrees of freedom 

Having discussed the physics of Eq. (1) in the absence 
of nuclear spins, we turn in this Section to the discus- 
sion of the effects of nuclear spin. One of the simplest 
Hamiltonians involving nuclear degrees of freedom would 
be realized in the case of one molecule per site: 



H = AY J SfT? 

i 

+i£> dd (R i -R i ) 



J %$4 S A 



-(SfS, 



As in Sec. II A 1, we ignored edge effects and dropped 
terms commuting with H. In this case, Tf is conserved 
on each site and becomes a classical variable. It can play 
the role of a tunable magnetic field at each site or the role 
of tunable disorder. The parameter A can be tuned to 
zero and away from zero, thus, decoupling nuclear spins 
from the rotor degree of freedom and coupling them. This 
tuning can be achieved, for example, by changing the 
strength of the DC electric field (see Sees. Ill and VA). 

While Eqs. (1) and (6) feature SfT? hyperfine inter- 
actions, other interactions between rotor and nuclear de- 
grees of freedom can also be generated. In particular, 
we show in Sec. VA that a judicious choice of rotor 
and nuclear states may allow for interactions of the form 
AS* T t z +A 2 Tl+A 3 Tf or even SfTr+S^Tf. Moreover, 
by combining the Hamiltonian in Eq. (1) or in Eq. (6) 
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with microwave and/or radio- frequency pulses applied at 
regular short intervals, one can use the Trotter approx- 
imation [77] to effectively modify the hyperfine interac- 
tion between Si and T$ from a simple SfT* interaction to 
any desired form. With this generalization, nuclear spin 
in Eq. (6) is in general no longer a classical variable. The 
Hamiltonian would then describe two types of spin-1/2 
species (each site having one of each): S species exhibit- 
ing interactions with neighboring sites and T species not 
exhibiting such interactions but interacting with S on 
the same site. Such a model is reminiscent of the Kondo 
lattice setup [78]. Moreover, the nuclear spin may allow 
to simulate the "orbital" degree of freedom, whose inter- 
play with spin (i.e. rotational) and charge (i.e. density) 
degrees of freedom may enable simulations of the exotic 
behavior of spin-incoherent Luttinger liquids [79] , transi- 
tion metal oxides [80] , and iron pnictide superconductors 
[81]. Finally, by using more than two nuclear spin states, 
one might be able, by analogy with alkaline-earth atoms 
[82, 83], to simulate exotic high-spin physics. 

While most of the discussion in the present manuscript 
focuses on quantum magnetism, the system also has 
promising quantum information applications. Having 
two outer electrons, alkali dimers have similar electronic 
structure to that of alkaline-earth atoms. Thus, one may 
consider extending some of the alkaline-earth quantum 
information processing proposals to polar alkali dimers. 
In particular, one can extend the idea of encoding quan- 
tum information in the nuclear spin degrees of free- 
dom from alkaline-earth atoms [84-90] to polar molecules 
[55, 91-93]. In this context, by analogy with alkaline- 
earth quantum register proposals [88] , information stored 
in the nuclear spins can be mapped via hyperfine inter- 
actions (or via microwave or radiofrequency fields) onto 
the rotor degree of freedom, which can then, in turn, be 
used to couple different molecules. By analogy with Ref. 
[88], we expect this system - particularly if more than 
two nuclear spin states are involved - to be useful in gen- 
erating high-fidelity many-body entangled states, such as 
cluster states or squeezed states. 



III. ROTATIONAL AND HYPERFINE 
STRUCTURE 

Having discussed in Sec. II the main features of Eq. (1), 
we present in Sees. Ill- VI the derivation of Eq. (1) and 
the ways, in which the coefficients in Eq. (1) can be con- 
trolled. Since we are interested in the effects of nuclear 
spin on the many-body Hamiltonian, we begin the deriva- 
tion of Eq. (1) by studying in this Section the rotational 
and hyperfine structure of a single molecule in the pres- 
ence of a DC electric field and zero or more continuous- 
wave (CW) microwave fields. The example molecules we 
are considering are 40 K 87 Rb [11] and 7 Li 133 Cs [15]. 

Following Refs. [17, 18, 20, 61, 94], the single-molecule 
Hamiltonian in the presence of a DC electric field and a 



CW microwave field is 

H = H + H mvj + Hhf, (7) 

where 

H = BN 2 - d E, (8) 
H = —r\ ■ ( K p p-iumwt ipp) 

H M = Hq + H IN + H t + H sc 

2 2 

= - e £ T 2 (VE 4 ) • T 2 (Q,) + CiN • li 

i=l i=l 

-c 3 y6T 2 (C)-T 2 (I 1 ,I 2 )+ C 4li-I 2 . (9) 

H describes the rigid rotor coupled to the DC elec- 
tric field. B is the rotational constant and N is the 
angular momentum operator describing the rotation of 
the molecule. The molecular quantization axis is chosen 
to be z, which is the direction of the applied DC elec- 
tric field (see Fig. 1). d is the dipolc moment operator, 
while d p = e p • d = dCp(9, (f>), where d is the permanent 
dipole moment of the molecule, p = 0, +1,— 1, and the 
spherical basis vectors are defined as eo = z and e±i = 

T-(x±zy)A/2 [61]. Here C*(6^) = ^aSl^M), 
where Yk tP are spherical harmonics, and spherical coor- 
dinates (6,<j)) describe the orientation of the rotor [61]. 

H mw describes the coupling of the rotor to a microwave 
field with amplitude E mw , frequency w mw , and polariza- 
tion e mw , which we assume to be equal to e_i, eo, or 
ei, which stand, respectively, for a~ , n, and a + polar- 
ization relative to the applied DC electric field. While 
H mw describes the action of a single microwave field, we 
will consider below the possibility of applying several mi- 
crowave fields, in which case i? mw would just feature the 
sum of the corresponding fields. 

Hm is the hyperfine interaction [17, 94], which is com- 
posed of four contributions: electric quadrupole Hq, 
spin- rotation Hxni tensor H tl and scalar H sc . These 
Hamiltonians couple the nuclear spins Ii and I2 of the 
two nuclei to N and to each other. The nuclei are num- 
bered as 1 = K and 2 = Rb for 40 K 87 Rb and 1 = Li and 
2 = Cs for 7 Li 133 Cs. 

The forms of H t and Hq warrant additional clarifica- 
tion. H t describes direct and indirect anisotropic inter- 
action between the two nuclei and is a scalar product 
[see Eq. (A2)] of two second-rank irreducible spherical 
tensors. The first tensor, T 2 (Ii,l2), is the second-rank 
tensor formed [see Eq. (Al)] out of Ii and I 2 . The second 
tensor, T 2 (C) = C 2 (#, </>), characterizes the orientation 
of the rotor and, hence, the relative position of the two 
nuclei. 

Hq describes the interaction between the electric 
quadrupole moment of each nucleus i and the electric 
field gradient at nucleus i due to the electrons and the 
other nucleus. Hq is also a scalar product of two second- 
rank irreducible spherical tensors. The first tensor is 
T 2 {Qi) = Qi ai^-i^ 'ft.Ii), where T 2 (T, T), is the 
second-rank tensor formed out of T and where eQi is the 



electric quadrupole moment of nucleus i. The second ten- 
sor is T 2 (VEi) = — ^-T 2 (C), where characterizes the 
negative of the electric field gradient at nucleus i. The 
values of all relevant molecular parameters for 40 K 87 Rb 
and 7 Li 133 Cs are given in Table I. All matrix elements 
are evaluated in Appendix A. 

At E = 0, the eigenstates of Ho are \N,M) obey- 
ing N 2 |A,M) = N(N + 1)\N,M) and N Z \N,M) = 
M\N,M). As we increase E, states with the same M 
mix to form the new eigenstates. Let us refer to the 
eigenstate that adiabatically connects to \N, M) (as we 
turn on E) as \4>n,m), as shown in Fig. 3(a). While 
\<Pn,m) are eigenstates of N z with eigenvalue M, they 
are not eigenstates of N 2 (for nonzero E); instead, they 
are superpositions of \N',M) for different N'. To al- 
low for a less cumbersome notation, let us also make the 
following simplifying definitions illustrated in Fig. 3(a): 
| A) = |(^jv,o) and \N) = |0jv,i)- In Sec. VB, we will also 
make use of the definition \N) = \<j)N,2) [see Fig. 3(a)]. 
The energies of \4>n.m) and the coefficients in the expan- 
sion of \<))n,m) in terms of \N',M) up to any desired N 
can easily be computed numerically by truncating the 
Hilbert space at some other - much larger - N. The 
fact that the splitting between \N, M) and \N + l,M) 
increases with N ensures that for any finite E, there will 
be some N above which the effect of E is negligible. 



A. Hyperfine structure in the simplest level 
configuration: {|mo), = {|0}, |1}} 

There is a great variety of possibilities - especially 
when microwave fields are applied - for choosing the two 
rotational states to play the role of |toq) and \mi) in Eq. 
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4U K 8V Rb 


v Li iM3 Cs 


d (Debye) 


0.566 [11] 


5.520 [95] 


B (GHz) 


1.114 [22] 


5.636 


B/d (kV/cm) 


3.9 


2.0 


d 2 / (47re (0.5Atm) 3 ) (kHz) 


0.39 


37 


h 


4 


3/2 


h 


3/2 


7/2 


{eQq) x (kHz) 


450 [22] 


18.5 


(eQq) 2 (kHz) 


-1410 [22] 


188 


ci (Hz) 


-24.1 


32 


c 2 (Hz) 


420.1 


3014 


c 3 (Hz) 


-48.2 


140 


c 4 (Hz) 


-2030.4 


1610 



TABLE I. Molecular parameters for 40 K 87 Rb and 7 Li 133 Cs. d 
is the permanent dipole moment, B is the rotational constant, 
and / is the nuclear spin. (eQq) characterizes Hq, c\ and C2 
characterize Him, C3 characterizes H t , and C4 characterizes 
H sc . In U, (eQq)i, and Ci = \p, the subscript i = 1 stands for 
K in KRb and for Li and LiCs, while the subscript i — 2 stands 
for Rb in KRb and for Cs in LiCs. The values for 40 K S7 Rb 
and 7 Li 133 Cs are taken from Refs. [17] and [20], respectively, 
unless otherwise indicated. 
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FIG. 3. (color online). (a) Eigenstates of Ho = 
_BN 2 — doE. (b-f) Level configurations employing mi- 
crowaves. The effective two-level dressed rotational degree 
of freedom {\m ),\m i)} is (b) {|0),VS|1) + V / T ^|2)}, (c) 
{|3),v^|l) + v / T^|2)}, (d) {^|0) + vT^|l),|l)}, (e) 
{v^|0) + v^T^ll),^!) + VT^b\(/>2,-i)}, (f) {Vo|2) + 
VT^a\2), Vb\l) + v/c|T> + v / l-fc-c|2)}. In figures (b-f), 
red (blue) levels make up the effective dressed rotor level |mo) 
(|mi». 



(1). In order to make the explanation of the main features 
of hyperfine structure clearer, we focus in this Section on 
the simplest example where no microwave fields are ap- 
plied and where |too) and \mi) correspond to the lowest 
two M — states: |mo) = |0) (= |</>o,o)) and |mi) = |1) 
(= 1 0i ,0) ) [see Fig. 3(a)]. In Sec. IIIB, we will extend 
this discussion to other level configurations. 

To simplify our effective Hamiltonian, we would like to 
prevent £/hf from coupling the states |0) and |1) to other 
states. Therefore, we need to assume that the applied 
DC field E is sufficiently large to split |1) from |1) and 
101. -1) by an amount larger than H^f. E.g., in KRb, to 
split off |1) from |T) and \4>i-i) by |(eQg)2|, one needs 
dE/B w 0.1. Since for KRb, B/d = 4 kV/cm, these 
values of dE/B are readily achievable. For LiCs, the 
required value of dE/B is even lower [dE/B = 0.015] 
since, for LiCs, \(eQq)2/ B\ is 40 times smaller. Moreover, 
in LiCs, B/d is 2 times smaller, which further reduces the 
required value of E. 

Under these assumptions, we can simply project H^f 
on the two states |0) and |1), without worrying about the 
crossterms: 



H 



m=0,l 



to) (to| (m|_ffhf \ m ) 



(10) 
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To understand the consequences of Eq. (10), let us 
follow the procedure similar to that in Ref. [20] and 
discuss what happens to different terms in H^f when 
we take the expectation value in a given rigid rotor 
state. (m\HjN\m) = for both states since these are 
both M = states and, therefore, give (m|N|m) = 0. 
{m\H sc \m) = H sc is unchanged since it does not involve 
rigid rotor coordinates. Considering Hq and H t , we have 
[using Eqs. (A14,A20)] 



(H Q ) = (P 2 (cos9))J2(eqQ) 



47,(2/, - 1) 



, (11) 



(H t ) = c 3 (P 2 (cos9)) U( 7 i%" + I1I2) - 2I?Jf)(12) 

Here P2(cos9) — Cq{6,4>) is the 2nd degree Legendre 
polynomial. Tg(Qi) acts on the i'th nucleus. We have 
used Eq. (Al) to get explicit expressions for T^I^Ij) 
and Tq(Ii,I 2 ). These expressions can also be obtained 
from Eqs. (A15,A21). 

Following Refs. [17, 21], we define the uncoupled basis, 
in which the two nuclear spin angular momenta are not 
coupled, and the coupled basis, in which they are coupled 
to form I = Ii +I2. The matrix elements are evaluated in 
Appendix A in both bases. We notice that H sc and (Ht) 
are diagonal in the coupled basis, while (Hq) is diagonal 
in the uncoupled basis. In Fig. 4, we plot (m\P2(cos9)\m) 
for m = 0,1,3 as a function of dE/B. An interesting 
"magic" point occurs at dE/B = 2.55: (O|P 2 (cos0)|O) = 
(l|P 2 (cos0)|l) = 0.18, i.e. the hyperfine structure in |0) 
and |1) is exactly the same. dE/B — 2.55 means 10 
kV/cm for KRb and 5 kV/cm for LiCs, so this point is not 
easy to access, but it could be useful for both quantum 
simulation and quantum computation applications as the 
point of decoupling of the nuclear and rotational degrees 
of freedom. As we can see in Fig. 4, similar "magic" 
points occur for the pairs of states {|0), |3)} and {|1), |3)} 
at dE/B < 10. 

Even if we are not at a "magic" point, where nuclear 
spin decouples from two rotor states |0) and |1), the hy- 
perfine structure is still relatively easy to understand. 
Hq competes with H sc to determine whether the un- 
coupled or the coupled basis is a good basis. Since in 



(m\P2(cos9)\m) 




dE/B 



FIG. 4. (color online). (m|P 2 (cos 0)\m) = (m\Cl (cos Q)\m) 
as a function of dE/B for m — (red), m = 1 (blue), and 
m = 3 (green). 



KRb (LiCs), (eQq)2 is 3 (2) orders of magnitude larger 
than C4, (rr^H^m) is almost diagonal in the uncou- 
pled basis, provided (m|P 2 (cos0)|m) > 10~ 3 (10~ 2 ). For 
example, the only place in the range of dE/B values 
shown in Fig. 4 where this condition breaks down for |0) 
(|1)) is near dE/B = 0(5), where (m|P 2 (cos0)|m) goes 
through zero. As pointed out in Ref. [20], at these points, 
(m\Hht\m) = -f^sc- Focusing for the moment on small 
values of dE/B, we find that in KRb (LiCs) H sc is dom- 
inant over H Q in |0) for dE/B < 0.1(0.5). On the other 
hand, we found above that we need dE/B > 0.1(0.015) 
to split |1) away from |1) by an amount larger than Hq. 
Thus, near dE/B = 0, there is a narrow range of dE/B 
for LiCs and no such range for KRb, where |1) is suf- 
ficiently split from |1), but H sc still dominates the |0) 
hyperfine structure. This observation supports the state- 
ment that for almost all values of dE/B, (rr^H^m) is 
dominated by Hq, while other hyperfine terms act as a 
perturbation. Therefore, in Figs. 5(a,b), we show the 
eigenvalues of ( Hq ) /(P 2 (cos 6)) for KRb and LiCs, re- 
spectively. From Eq. (11), we see that these eigenvalues 

are Yli=i ( e lQ)i ^^I'j. (2/ - ' wnere is the magnetic 
quantum number of nucleus i. 

H sc and H t can then be treated as a perturbation: 



(H IN ) = (c 4 - 2c 3 (P 2 (cos(?)))/ 1 % 2 
(P 2 (cos0)))(l+l 2 -+lrl+). (13) 



H sc + (H t ) 
V 2 2 

Here If If is diagonal in the uncoupled basis and just 
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FIG. 5. (color online). Eigenvalues (in kHz) of 
(Hq)/(P 2 (cos0)) for KRb (a) and LiCs (b). The horizon- 
tal axis is the magnetic quantum number Mi (for the K nu- 
cleus in KRb and for the Li nucleus in LiCs), while \M 2 \ is 
indicated separately for each group of levels. The two red 
circles indicate (Mi,\M 2 \) = (0,1/2) and (1,1/2) in (a) and 
(Mi, |M 2 j) = (-1/2, 1/2) and (1/2, 1/2) in (b). 
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shifts the energies slightly. The flip-flop term ij 7 2 + 
IfiJ changes (Mi,M 2 ) by (1,-1) or by (-1,1). This 
term is important provided the energy difference between 
the two states it connects is smaller than the flip-flop 
strength. For a typical value of (P^cos^)) ~ 0.1 (see 
Fig. 4), the smallest relevant splitting in (Hq) for KRb 
is between {M x , M 2 ) = (0, 1/2) and (1,-1/2) [red circles 
in Fig. 5(a)] and is equal to about 1 kHz. Since C4 in 
KRb is « —2 kHz, a few of the states in Fig. 5(a) will get 
mixed by H sc + (H t ), but for most states the uncoupled 
basis stays a good basis. The situation is similar in LiCs, 
where the smallest relevant nonzero splitting in (Hq) is 
between (1/2,1/2) and (3/2,-1/2). For (P 2 (cos<9)) - 0.1, 
this splitting is equal to ~ 1 kHz, which is comparable 
to C4 = 1.6 kHz. An additional feature in LiCs is that 
(Hq) has degenerate states (-1/2,1/2) and (1/2,-1/2) [red 
circles in Fig. 5(b)] that get mixed by the flip-flop term. 

It is worth pointing out that the application of a mag- 
netic field can help in defining the uncoupled basis as a 
good basis, as is done in current experiments [22]. For 
example, this knob can be used to make the above dis- 
cussed flip-flops off-resonant or to decouple the nuclear 
spins when (P 2 (cos#)) is small. In particular, this allows 
one to decouple the nuclear spins from each other in state 
|0) even at small DC electric fields [22]. 



B. Hyperfine structure in other level 
configurations 

In the previous Section [Sec. Ill A], we described the 
hyperfine structure in the simplest configuration of rota- 
tional levels: {|mo),|mi)} = {|0),|1)}. In this Section, 
we extend this discussion to other configurations of rota- 
tional levels. 

Even without microwave fields, a great variety of possi- 
bilities exist for choosing rotational states to prepare the 
effective rotor degree of freedom that is featured in Eq. 
(1). One could, for example, choose, instead of states 
{|0),|1)}, the states {|1),|3)}, which, as we will see in 
Sec. VI, may have some advantages over the former. 

One could also consider {|rno),| m i)} = {|0),|1)} or 
{|mo), l^i)} = {|1), |1)} as the effective rotor degree of 
freedom. In order to avoid the coupling of |1) to |</>i,_i) 
by Hq and H t (and later by the dipole-dipole interac- 
tion), we can apply, for example, a er~-polarized mi- 
crowave field coupling \4>i,-i) to |0 2 ,— 2) that would shift 
the state . Once this is done, H^f can be projected 

on each of the two states as in Eq. (10). Another impor- 
tant difference will be the fact that (1|N|1) = z 7^ 0, so 
that (l|-ffjjv|l) = J2i c i^i- This term will contribute to 
(l|iJhf|l) in Eq. (13). Being diagonal in the uncoupled 
basis, the term (1\Hjn\1) will just slightly shift the lev- 
els obtained after diagonalizing (1\Hq\1). This term may 
provide an extra control knob. In particular, in LiCs, c 2 
is about twice the value of the scalar coupling C4 and 
will, thus, play an important role for nuclear spin states 
that are nearly degenerate under (Hq) . In addition to 



being a control knob, Hjjq may also give rise to some 
complications. Specifically, the point where (1|P2|1) is 
equal to ( 1 1 J-2 1 1 ) is, in fact, not an exact magic point for 
the two states (i.e. the two hyperfine structures do not 
perfectly match) due to the Hjn term. However, first, Cj 
are rather small (a few orders of magnitude smaller than 
the dominant quadrupolar term - see Table I). Second, 
( 1 1 i^iN 1 1) can vanish exactly for If = If = 0, which can- 
not happen for our isotopes but is, in general, possible. 
Third, one can slightly adjust the value of the DC electric 
field from the one that gives (1|P 2 |1) = (IIP2II) m such 
a way that some (but not all) desired nuclear spin states 
have the same relative energies in |1) and |1). 

The application of microwave fields allows to gain bet- 
ter control over the effective Hamiltonian [21, 26-28, 30- 
32, 36, 48, 61-64]. In this Section, we consider two exam- 
ples of microwave control. In the first example, proposed 
in Ref. [32], we couple states |1) and |2) with a linearly 
polarized microwave [Fig. 3(b)]. We assume that the mi- 
crowave field is sufficiently weak that it can be treated 
within the rotating-wave approximation and that its off- 
resonant couplings on other transitions can be ignored. 
Furthermore, we assume that the microwave Rabi fre- 
quency f2 = -E mw (2|G?o|l) is much larger than the hyper- 
fine structure splittings, so that all hyperfine transitions 
are addressed equally. In principle, weak microwave fields 
coupling individual nuclear spin levels can also be used to 
implement quantum magnetism with polar alkali dimers 
[22, 26, 27]; however, for simplicity, we will not discuss 
this case in the present manuscript. In KRb, a Rabi 
frequency spanning all hyperfine levels (of order a few 
MHz) requires a microwave intensity of a few W/cm 2 , 
which is achievable in the laboratory. In LiCs, which has 
a 10 times larger dipole moment and 10 times smaller 
hyperfine splittings, the required microwave intensity is 
10 4 times smaller. The two requirements of staying off- 
resonant with other rotor transitions (and staying within 
the rotating-wave approximation) but at the same time 
addressing all hyperfine levels can easily be achieved since 
in KRb (LiCs) the splitting between the rotor levels ~ B 
is 3 (4) orders of magnitudes larger than the largest hy- 
perfine constant (eQq)2 (see Table I). 

The application of the microwave field will produce, in 
the rotating frame, two dressed states [63] . One of them 
will form the state \m\) = \/a\l) + \A — a|2) , where we 
assumed for simplicity real positive coefficients and where 
a can be controlled by the amplitude and detuning of the 
microwave field. Projecting Hm on the subspace spanned 
by \mo) = |0) and |mi), we obtain the following form of 
the hyperfine interaction 

flhf« |0)(0|(0|ff h f|0) + 

+|mi)<mi|(a(l|fl hf |l) + (1 - a)(2|P hf |2)).(14) 

Notice that (l|Phf|2) does not contribute since, in our 
rotating frame, it is rapidly oscillating. The discussion 
of the {|m ),|mi)} = {|0),|1)} configuration then ap- 
plies with the change that (l|Phf|l) is replaced with 
a(l|Hhf|l) + (1 - a)(2|f/hf|2). One advantage of this 
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configuration over the {|0), |1)} configuration is that, for 
a given choice of the DC electric field, the magic point 
where (0|P 2 |0) = a(l\P 2 \l) + (1 - a)(2|P 2 |2) may be ac- 
cessed by tuning a. At this magic point, the nuclear and 
rotational degrees decouple, as discussed above. 

The second example involving microwave fields that we 
consider in this Section involves the application of a a + 
microwave field near resonance with the |0) — |1) tran- 
sition. We pick one of the two rotating-frame dressed 
states |mo) = \/a\0) + \f\ — o|l) as one of the two effec- 
tive rotor states and state \m\) — |1) as the other [Fig. 
3(d)]. In contrast to the microwave-free {|0), |1)} config- 
uration, in this example, we can safely ignore the state 
|</>i, -1) assuming the dressed state \mo) is shifted by the 
applied microwave sufficiently far away from the state 
|0i ( _i). If \mo) is too close in energy to \(pi-i), then 
|0i -l) can be shifted away using a separate microwave 
field coupling it, for example, to |02,- 2 ). Projecting the 
hyperfine Hamiltonian onto |mo) and |1), we obtain the 
same Hamiltonian as in Eq. (14) except states |0), |1), 
and 1 2) get replaced with states |1), |0), and |1), respec- 
tively. 

Of course, numerous other coupling schemes are also 
possible. For example, one can apply two microwave 
fields acting on two different transitions and use one 
dressed state from each transition as the basis. One can 
even consider appling more microwave fields, as we will 
do in Sec. VB. The above discussion of the hyperfine 
structure can be readily extended to these cases. 



IV. OPTICAL POTENTIAL AND TENSOR 
SHIFTS 

As discussed in Refs. [21, 27, 30, 61, 96], a rigid rotor 
placed into an optical lattice experiences level shifts - 
called tensor shifts - that depend on the internal state of 
the rotor. In this Section, we summarize the derivation 
of tensor shifts from Ref. [61], consider ways to control 
these shifts, and discuss the effects of these shifts on our 
Hamiltonian. 

Following Ref. [61], we consider an off-resonant light 
field E opt (R,t) = E opt (R)e- iw * + ex.. We recall that 
we use the X-Y-Z coordinates to describe the 2D optical 
lattice, which lies in the X-Y plane, while the x-y-z co- 
ordinates will have z along the applied DC electric field 
(see Fig. 1). In the present Section, we will ignore the 
hyperfine structure - we will put together the optical 
potential and the hyperfine structure in Sec. V. The AC 
Stark shift Hamiltonian acting on a rigid rotor describing 
the ground electronic and vibrational state of a molecule 
is then 



flopt(R) = -Eopt(R)* • &(u) ■ Eopt(R), (15) 



where 



a((jj) — a±(oj) (16) 
+ [anH - ax(w)] £(~l) p C_ p (#, 0)^,(0, 0)e p ® e* pl . 



Here (9,4>) are the spherical coordinates of the rotor. 
a||(w) and a±(uj) are dynamical polarizabilities at fre- 
quency uj parallel and perpendicular to the rotor axis. 
Due to the difference in matrix elements and in the en- 
ergy difference between states contributing to the two 
polarizabilities, an and a± are generally different giving 
rise to the term oc [0:11(0;) — a_i_(a;)] describing a rotor- 
state-dependent shift [second line in Eq. (16)]. 

We suppose that E opt (R) = E(R) J2 p =~i Pp& P , where 

S P Pp&p i s a un it vector (i.e. J2 P \fip\ 2 = 1) describing the 
polarization of the light, which, for simplicity, we assume 
to be spatially uniform. We then find 



ff p t (R) = -|ff(R)| 2 [ M > 3 111 ' 
2 

+ [a N (w) - a±(u)] 7pCp(6»,(5 



(17) 



P^i(3±i, 7±i 



v/3 



where 7± 2 = - 

7o = |/3o| 2 -|. 

We now recall that we will be working at DC electric 
fields that are large enough to separate the rotor states 
of interest from all the other states by a shift larger than 
the hyperfine interaction strength (> 1 MHz). In the 
cases where the state |1) is involved, we assume that a 
microwave field acting on |1) itself or on 101,-1) splits 
the two by a similarly large shift. Since 1 MHz is greater 
than typical optical lattice potential strength (10 — 100 
kHz), the lattice potential is too weak to induce tran- 
sitions between the rotor levels, and we can therefore 
just evaluate Popt in each rotor state. Moreover, any 
l m ) = |0JV,m) (with any M) is an eigenstate of N z , so, 
for p^O, (m\Cp(8, 0)|m) = 0. Therefore, for such states 
|m), we get the microwave-free optical potential 



ff opt (R) = -|£(R)| 2 [a H- 
+o 2 (u;) y^(m|P 2 (cos 9)\m) \m) (m\ 



(18) 



where 



oc (uj) 



2a±(uj) + 0|| (uj) 



a 2 (w) = [a\\(uj) - a ± (u)} (\(3 \ 2 - i 



(19) 



p,p' 



The dependence of tensor polarizability a 2 (u>) on f3 is 
in direct analogy with the corresponding dependence in 
atomic tensor polarizabilities [97, 98]. 

In the case where a microwave field is applied, the 
optical potential can be computed as follows. For the 
{\m ), |mi)} = {\/a|0) + y/l - a\l), |1)} configuration 
[Fig. 3(d)], the optical lattice potential is 

ff op t(R) = -|^i (R)| 2 [«oM + o 2 (u; 1 )(|l)(l|(l|P 2 |l) 
+ |m o )(m o |(a(0|P 2 |0) + (1 - a)(T|P 2 |T)))l . (20) 
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Up to a relabeling of states, a similar expression holds 
for the {|0), y/a\l) + VT- ^a\2}} configuration [Fig. 3(b)]. 

By changing the frequency and the polarization of the 
applied light and by using several [27, 61] laser beams at 
different frequencies or polarizations, one can control the 
strength of the tensor shift relative to the scalar shift. In 
particular, it is often convenient to have a lattice that 
is independent of the rotor state. One can envision the 
following avenues for achieving this. 

First, as already pointed out in Ref. [96], for any pair of 
states m and m', the tensor shift vanishes at the "magic" 
points in Fig. 4, where (m\P2(6)\m) = (m'\P2(9)\m') [see 
Eq. (18)]. We recall that these are the same points where 
the nuclear spins and the rotor degree of freedom decou- 
ple. In the case where a microwave field is applied, one 
has an extra control knob to arrive at the "magic" point 
for the two states of interest. For the example consid- 
ered in Eq. (20), the microwave Rabi frequency and de- 
tuning can be used to control a to obtain a lattice that 
is the same for states |mo) and |1), which happens when 
a{0|P 2 |0) + (l-a)(T|P 2 |I) = (l|P 2 |l). 

Second, by analogy with "magic" frequencies for 
atomic levels [97, 98] and for vibrational molecular lev- 
els [99], one may look for a "magic" frequency oj, for 
which c\!||(w) = a±(uj), in which case a 2 (w) would van- 
ish. However, the search for such a "magic" frequency 
may be significantly complicated by the requirement to 
keep spontaneous emission low [100]. 

Third, as already pointed out in Ref. [96], a 2 (uj) would 
also vanish if one chooses a polarization, such that \/3q\ 2 = 
1/3. For example, a linear polarization making an angle 
cos _1 (l/v3) with the z-axis (i.e. with the DC electric 
field) would work. 

Fourth, one may use two laser beams [27, 61] that 
have a 2 of opposite signs. Assuming these beams can be 
made to have the same spatial profile (which can be done, 
for example, with holographic techniques [101] or angled 
beams [102]), their relative intensities can be adjusted in 
such a way that the combined tensor shift vanishes. 

Finally, if a 2 of opposite sign is difficult to achieve, 
as long as a 2 /ao is different for the two lasers, one can 
choose the two lasers (on the example of ID) to have 
spatial profiles £ 2 cos 2 (XX) and Xf sin 2 (XX), respec- 
tively (for some wavevector X). By tuning the relative 
intensities of the two lasers, one can achieve X 2 a 2 (wi) = 
E\a.2^2) (where ui\ and w 2 are the laser frequencies), 
which would allow to make the tensor shift spatially in- 
dependent [cos 2 (XX) + sin 2 (XX) = 1]. The spatially 
independent shift can then be treated as a slight modifi- 
cation to the internal structure. While this last solution 
described a ID lattice, three ID lattices can be combined 
into a 3D lattice provided their frequencies differ slightly, 
so that the lattices do not interfere. 



V. DERIVATION OF THE HAMILTONIAN 

In this Section, we use the results of Sees. Ill and IV to 
derive the Hamiltonian in Eq. (1) and to show how var- 
ious terms in this Hamiltonian can be tuned. We recall 
that, as shown in Fig. 1, the molecules are confined to the 
X-Y plane and are subject to a 2D optical lattice in that 
plane. We also recall that a DC electric field of strength 
E is applied in the direction z that makes a polar angle 
0o with the Z-axis and has an azimuthal angle $ m the 
X-Y plane. The system is then described by five one- 
body Hamiltonians and one two-body Hamiltonian. The 
five one-body Hamiltonians are [17, 18, 20, 61, 94] 



H Q = BN 2 - d E, 
H^mw — d • (X mw e mw e -\- c.c.^ 

Hm = Hq + Hin + H t + H sc , 
H opt = -E p t (R)* • a(cu) • E pt(R), 



P 

2Mn 



(21) 
(22) 
(23) 
(24) 

(25) 



The molecules are assumed to be in the electronic and 
vibrational ground state, Hkin describes the kinetic en- 
ergy, and M m is the mass of the molecule (subscript m 
here stands for the word molecule to avoid confusion with 
the magnetic quantum numbers M). 

The two-body Hamiltonian for molecules 1 and 2 is 
given by the dipole-dipole interaction 



H&d = 



1 



4vrenP 3 



d«-d( 2 )-3(R-d«)(R -01^)1(26) 



Here d( J ) is the dipole moment of molecule j and R = 
HR, is the vector connecting the two molecules. 

From Table I, we see that the KRb system has a con- 
venient separation of energy scales, which we have al- 
ready used in Sees. Ill and IV: Hq ~ B ~ 1 GHz, 
H hi ~ H Q ~ 500 kHz, H opt + #ki„ - 10 - 100 kHz, 
Hdd ~ d 2 / (4ireoR 3 ) ~ 1 kHz (where we assume a typi- 
cal separation R ~ 0.5//m between two neighboring sites 
of an optical lattice). As discussed in Sec. Ill, we also 
choose H mw to have an energy scale significantly below 
Hq and significantly above H^f. Therefore, the Hamilto- 
nians can be treated in order of decreasing energy scale. 

In the case of LiCs, dipole-dipole interactions [see Ta- 
ble I] are typically on the same order or even stronger 
than the optical potential. In that case, the physics 
changes and involves such effects as Wigner crystalliza- 
tion [55, 62]. Wigner crytallization has the exciting po- 
tential of bringing the molecules closer together (for ex- 
ample, if the optical lattice is not present) and producing 
strong internal-state-dependent interactions. However, 
the study of such models involves the phonon modes [55] 
and is beyond the scope of the present work. Therefore, 
in the present Section, we assume that we either work 
with KRb or that the rotational levels of LiCs are cho- 
sen in such a way [for example, using states |1) and |3) 
at small DC fields (see Fig. 7) or employing microwaves] 
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that dipole-dipole interactions are much weaker than the 
optical potential. 



A. Derivation of the Hamiltonian for the simplest 
level configuration: {|mo}, = {|0}, |1}} 

To derive the Hamiltonian in Eq. (1), let us begin in 
this Section with the simplest case where no microwave 
fields are applied and where we restrict ourselves to rotor 
states |m ) = |0) (= |<Ao,o)) an d |mi) = |1) (= |</>i,o))- 
We will consider other level configurations in Sec. VB. 

The diagonalization of H + Hm + i? mw was discussed 
in Sec. III. In particular, we showed that for a generic 
DC electric field, the hyperfine structure in |0) and |1) is 
almost diagonal in the uncoupled basis. Therefore, if we 
would like to ignore the nuclear spin, one way to do this 
is to prepare all molecules in a nuclear spin state that is 
an eigenstate of both (0|iJ M |0) and (l\H M \l). The hy- 
perfine energy can slightly change the energy difference 
between |0) and |1), but the total number of molecules in 
|0) and the total number of molecules in |1) will be sepa- 
rately conserved, making the precise value of the energy 
difference between |0) and |1) unimportant. Another way 
to ignore the nuclear spin is to use the "magic" point at 
which the hyperfine structure of |0) and |1) is exactly 
the same, in which case one does not even have to pre- 
pare all molecules in the same nuclear state to observe 
nuclear-spin-independent dynamics. 

To include the nuclear spin into our dynamics in a min- 
imal way, we pick two nuclear spin states | f) and | \.) 
that are eigenstates of both (0|£/hf|0) and (l|iJhf|l)- This 
can easily be done since, for a generic DC electric field, 
the two hyperfine structure Hamiltonians are almost di- 
agonal in the same (uncoupled) basis [see Sec. III]. In the 
second quantized notation, we can, thus, write 



hf 



ma 



rll. 



(27) 



where n mtT is the number of molecules in internal state 
m(=0,l), a(=t,|). 

We now consider H opt + H^ n . We suppose that the 
molecules are confined to the lowest band of a 2D lat- 
tice in the X-Y plane with the third direction Z frozen 
out. As discussed in Sec. IV, |0) and |1) will generi- 
cally feel lattices of different strength, so that H opt = 
Yl m =o i l TO )( TO |Kn(R)- We can then expand the molec- 
ular operator \E'„ lcr (R) in (real) Wannier functions as 
\l/ mcr (R) = J2j Wj m (R)cj mc7 , where j sums over sites in 
the X-Y plane. Here Wj m (R) = w m (R — Rj), where 
Rj is the position of site j in the 2D lattice. Absorbing 
zero-point energy into Eq. (27), -ff pt +-ffkin can then be 
rewritten as 



H op t + -ffkin 



(ij)ma 



C lmcr C jmc + n - C - 



(28) 



where the sum is taken over all nearest neigh- 

bor pairs and where the tunneling amplitudes 



- / d 3 -Rw. lm (R)[-V 2 /{2M rn ) + V m (R)}w jm (R) for i and 
j nearest neighbors. For simplicity, we assumed that tun- 
neling amplitudes are the same for all nearest neighbor 
pairs. 

We now consider H dd . Since both |0) and |1) are M = 
states, (extended to many molecules) simplifies to 
[see Eqs. (A25-A27)] 



(29) 



In second quantized notation, and keeping only the 
terms that conserve the total number of molecules in 
state m (for each m), H^d can be rewritten as 



H dd = \Y,J d3Rd3R "MR - R') 



(30) 



-{^-^'*m CT (R-)^V'(R')*mV'(R')*m CT (R) 



-Mo 



x [*S <r (R)*L(B- , )*Oa'(R , )*la(R)+h.C.] } 



where fJ, m m' = ( m l^o| m ') is the transition dipole moment 
between \m) and \m') and where fi m — (m\do\m) is the 
dipole moment of state |m). The presence of nonzero 
dipole moments no and fj,\ is expected since |0) and |1) 
are eigenstates of the rigid rotor Hamiltonian in the pres- 
ence of a DC electric field. One should keep in mind that 
certain values of dE/B may give rise to terms that do not 
conserve the total number of molecules in state \m): for 
example, at dE/B k, 3.24, the energy difference between 
|0) and |1) is equal to the energy difference between |1) 
and |2) = |</>2,o), and dipole-dipole interactions can reso- 
nantly turn two molecules in state |1) into a molecule in 
state |0) and a molecule in state |2). We assume, how- 
ever, that we avoid such accidental degeneracies. 

Expanding 5 , mo .(R) in Wannier functions, we obtain 

H dd = \ E / rf 3 Rrf 3 R'Vdd(R - R') 

hnnu 

mm' a a 1 

m^m'Cj ima Cj 2rn , a ,Cj 3m ' a 'Cj 4ma 

l - J d 3 Rd 3 R'V d d(R-R')w n o(R)w j2l (R') 



3U23334 



cw J 3 (R')w 34 i(R)Moi c ]iO<t c ]21<t' c ^ 0ct ' c ^ 1ct + h ' c - 



(31) 



Here the hardcore constraint means that j\ ^ ji an( i 
3s 7^ J4- We now make two approximations: (1) the 
extent of w is much smaller than the distance between 
the sites, and (2) only terms where i = j\ = 34 3 ; = 
32 = 33 contribute. These approximations allow to take 
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Vdd(R-i ^ Rj) outside of the integral. The result is 
H dd = 9 X! ^dd(Ri - R-j) 



E 



j± 



, (32) 



where Jj_ = 2/Xqi (/ d 3 RiUio(R)wii(R)) 2 - Interestingly, 
the presence of tensor shifts, thus, does not affect the co- 
efficients of n im rij m > because the Wannier functions are 
always normalized. The only effect of tensor shifts is, 
thus, a slight reduction of Jj_ from its tensor-shift-frce 
value of The latter makes perfect intuitive sense: 

the matrix element is reduced due to reduced overlap. 
One can view this effective modification of /^oi as an ex- 
tra control knob. In the remainder of Sec. V, however, 
we will assume for simplicity that all rotor states feel 
the same optical potential; tensor shifts can easily be in- 
cluded by analogy with the above example and will lead 
to similarly reduced matrix elements. Expressing rii m 
in terms of rij and Sf , we find that Eq. (32) is equiva- 



-, w 



lent to -ffdd in Eq. (1) with V — , . ,. — 
and J z = (/Lto — A*i) 2 - In Appendix B, we calculate cor- 
rections to the approximations made to arrive at Eq. 
(32). While these corrections lead to interesting effects, 
such as interaction-assisted tunneling, these corrections 
are small. It is worth pointing out that at an electric 
field of dE/B = 0.1, which we need to prevent i?hf from 
coupling |1) to |1) and \<fri-i) and to decouple the nu- 
clei in state |0), Moi/^o ~ 300. This means that at this 
value of dE/B, the values of V, W, and J z are negli- 
gible compared to J±, making the V = W = J z — 
model studied in Ref. [7] applicable. In row #1 of Table 
II, we collect the values of V, W, J z , and J± (for the 
case of no tensor shifts) and list the main features of the 
{\m ), |mx)} = {|0), |1)} scheme. 

Let us now simplify the internal state Hamiltonian in 
Eq. (27). Using the definition n ima = cl ma c ima , Eq. (27) 
can be rewritten as 



H + H : 



hf 



E. 



ma ' L imo • 



(33) 



We will use conservation laws to simplify this expression. 
In particular, our Hamiltonian conserves the total num- 
ber of molecules (no = not + n oi)i the total number 
of 1 molecules (ni = ni-j- + n±i), the total number of 
j" molecules (n-j- = not + nif), and the total number of 
1 molecules (n± = n t + ^14,)- Only three out of these 
four quantities are independent since the first two and 
the last two quantities both sum to the total number of 
molecules. Thus, subtracting from the final Hamiltonian 
constant quantities that commute with it, the only rele- 
vant internal-state Hamiltonian will be 

Ho + H hi ->• -(riiot - riioi - n in + rim) 



i 



(34) 



where we assumed that there is at most one molecule per 
site. Here 



A = 



((O|P 2 (cos0)|O)-(l|P 2 (c 



(35) 



where the last approximation is made provided Hq dom- 
inates the hyperfine structure and where | f) = \Mi, M 2 ) 
and I I) = \M[,M' 2 ). We see thus that A can be tuned 
with a significant degree of flexibility. In particular, Fig. 
4 shows that (O|P 2 (cos0)|O)-(l|P 2 (cos0)|l) can be tuned 
by adjusting dE/B. On the other hand, Fig. 5(a) shows 

(on the example of KRb) that J^Li 4j/(2/?-i) K M ») 2 " 
(Ml) 2 } can be adjusted between 12 kHz [e.g. (M[,M 2 ) = 
(0, 1/2) and (M U M 2 ) = (1, 1/2)] and - 1 MHz. 

Let us now briefly discuss the possibility of obtaining 
more complicated interaction terms between Si and T; 
than the simple AS*T? in Eq. (34). First, it is possible 
to get a Hamiltonian of the form ASfT? + A 2 T? + A 3 T? '. 
In the case of one molecule per site in the absence of 
tunneling, such a Hamiltonian still conserves Sf as the 
original ASfT* Hamiltonian but no longer conserves Tf. 
The term Tf can be obtained by working in the regime 
when the term I 2 + If I 2 in Eq. (13) couples the two 
chosen spin states and is not negligible. Whenever Tf 
is not negligible, the term A 2 T? arrises naturally fol- 
lowing a derivation similar to that leading to Eq. (34). 
Second, it is also possible to obtain terms of the form 
S+T[~ + 5 4 ~ Tj + . Such terms allow one to exchange S and 
T excitations within the same molecule. We can obtain 
such terms by using, for example, states |1) and |0i,-i) 
as the two rotor states. In that case Hq and H t can 
cause transitions between these two levels while at the 
same time changing 7f + If by 2. 



B. Derivation of the Hamiltonian for other level 
configur at ions 

In the previous Section [Sec. VA], we derived the 
Hamiltonian in Eq. (1) for the simplest level configu- 
ration: {|too), = {|0),|1)}. In this Section, we 
show how the coefficients V, W, J 2 , and Jl in Eq. (1) 
can be controlled by choosing other level configurations, 
including those configurations that involve one or more 
microwave fields. 

The results are summarized in Table II. The 
microwave-free {|1), |3)} scheme (#2 in Table II) has 
the same form of the dipole-dipole coefficients as the 
{|0),|1)} scheme. However, it has two important fea- 
tures that distinguish it from the {|0), |1)} scheme: the 
transition dipole moment between |1) and |3) vanishes 
for E = 0, and the permanent dipole moments of states 
|1) and |3) point in the same direction at small fields E. 
This may help stabilize the system against chemical re- 
actions (see Sec. VI) and may help reduce the strength 
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R,otor st cites used 


Expressions for V, W 7 , J z , Jx 

V r =po + Ai) 2 +So + -Bi]/4 
W = [Ao + So — A 2 — Bi]/2 

Jz — (AO — sii ) ■+■ £)() ■+■ £>i 


Special features 


#1 


|m ) = |0> 

K) = |i> 

Fig. 3(a) 


y = ( M0 + Ml )74 

W = (Mo - M?)/2 
J z = (mo ^ Mi) 2 
Jx = 2Moi 


- Simplest. 

- At small dE/B, V « W » J z « and J ± > 0, 
yielding the dipolar t-J± Hamiltonain [7]. 


#2 


\mo) = |1) 
jmi) = |3) 
Fig. 3(a) 


V = (mi+M3)V4 
W={j&- Ml)/2 
Jz = (mi - M3) 2 

T r» 2 

Jx = 2m1 3 


- At small dE/B, MiM3 > Ml3i which may help 
stabilize the system against chemical reactions 
(see Sec. VI). 


#3 


\mo) = |0> 
|mi) = |1) 
Fig. 3(a) 


1/ = (mo + Mt)V4 
W = (Mo - Mt)/2 
Jz = (mo - mt) 2 

T 2 

J ± = -Mot 


- Simplest configuration with J± < 0. 

- A microwave field is required to shift \4>i,-i) 
out of resonance with 1). 


#4 


jm > = |0) 

jmi) = Va|l) + Vl - a\2) 
Fig. 3(b) 


Ao = Mo 

Ai = aMi + (1 — a)M2 
Bo = 

Si = 2M? 2 a(l - a) 

Jx = 2(Moia + M02(l - a)) 


- At (dE/B, a) = (1.25, 0.74), W = 0, 
Jz = Jx = 0.36d 2 , and V = O.lJz, making 
Eq. (5) similar to the SU(2)-symmetric 
t-J-V model, which exhibits suppressed phase 
separation [58]. 


#5 


|m ) = y/a\0) + VI ~ o|X) 
|mi) = |1> 
Fig. 3(d) 


A = aMo + (1 — a)mt 
Ai = mi 

-Bo = -MoT a ( 1 - a ) 
Si = 

Jx = 2<jmoi - (1 - a) MiT 


- Jx = can be achieved at any dE/B by 
adjusting a. 

- J z < can be achieved. 

- V < can be achieved. 


#6 


|m ) = |3> 

|mi) = a /o|1)+a/1-o|2) 
Fig. 3(c) 


A = M3 

Ai = a/xi + (1 — a)M2 

5 = 

51 = — a(l — ffl)Mi3 

Jx = 2aM? 3 - (1 - a)M 2 2 


- J z — and Jx = lines intersect in (dE/B, a) 
space at (dE/B, a) = (2.6,0.92). So if 
we write J z = \J\ cos if) and Jj_ = | J| sin-i/i, 
arbitrary ip can be achieved around that point. 


#7 


jm ) = a|0) + VI - a\l) 
|mi) = 6|1> + Vl — 6]02,_i) 
Fig. 3(e) 


A = omo + (1 — ci)mt 
Ai = 6/xi + (1 — b)M 2 

5 = -a(l - o)Mot 

51 = -6(1 - 6)m?3 

Jx = 2a6Moi - a(l - 6)mq 5 
-(1 -«)&Mti 


- At (dE/B, a, 6) = (1.7, 0.33, 0.81), W = and 
Jz = Jx = -4V = 0.089d 2 , making Eq. (5) 
very similar to the standard t-J model [8]. 


#8 


|m > = Vo|2) + VI - a\2) 

mi) =V6|1) + Vc|l> 

+Vl-6-c|2) 
Fig. 3(f) 


Ao = dM 2 + (1 — ci)m 2 

Ai = 6mi + cmt + (1 — 6 — c)m2 

5 = -a(l - a)M|j 

51 = -bcfi 2 lT - c(l - 6 - c)m2t 

+26(1 - 6 - c)m? 2 
Jx = 2(1 - a)cM§j - (1 - a)6M§j 

-(1 — a)(l — 6 — c)m§ 2 — ac M^y 


- The manifolds V = 0, W = 0, J z = 0, 
and Jx = intersect in (dE / B, a,b, c) space 
at (dE/B,a,b,c) = (2.97,0.059,0.56,0.38). 
Full control over V, W, J z , and J± is 
achievable around that point. 



TABLE II. The expressions for the dipole-dipole interaction coefficients for several different level configurations. For config- 
urations #1 through #3, the coefficients V, W, J z , and J± are listed directly. For other configurations, we instead list the 
expressions for Ao, Ai, So, Si, and J±; the expressions for V, W, and J z can be computed from A v and S p using the formulas 
provided at the top of the table. While the presented expressions for the interaction coefficients assume no tensor shifts, the 
effect of tensor shifts is straightforward to include. Some notable features of each configuration are noted in the last column, 
while a more detailed discussion is provided in the text. 

of dipole-dipole interactions in LiCs below the strength 
of the optical lattice potential, which is necessary for the 
applicability to LiCs of the treatment that we present. 

To calculate the coefficients V, W , J z , and J± in the 
{|0), |1)} scheme (#3 in Table II), we have to extend Eq. 
(29) to account for the fact that d + and d- now play a 
role (recall that d± = e±i -d). However, energy conserva- 
tion still forces the conservation of total N z of the two in- 



teracting molecules, thus, making sure thatT2(dW,d«) 
contributes only for p — [see Eq. (A25)]. Therefore, 
according to Eq. (A26), d^d^ in Eq. (29) should be re- 
placed with 4 i) 4 j) + U d + d( -+ d - )d + > )- Projecting on 
the states |0) and |1) and ignoring off-resonant terms, we 
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obtain 



(|0} i (0| iM o + |l)i<l|iMT) (|0>i<0| J A* + |1>j<1|jMt) 



|01)(10| +h.c. 



(36) 



Here fj, m — (m\do\m), fi Q j — (0|d_|l), and \mm!) means 
that molecule i (j) is in state \m) (|m')). From now on, 
we will use the natural notation that /i mm < — {m\d p \m!) 
is the transition dipole moment between \m) and \mf) 
computed using d p for the appropriate p. The interac- 
tion in Eq. (36) has the same form as the corresponding 
interaction for the {|0), |1)} scheme except /z^ is replaced 
with —fi^j/2. This minus sign comes from the physical 
effect that two dipoles rotating in the x-y plane give an 
averaged interaction that is equal to negative one-half of 
the interaction for two dipoles pointing in the z direction 
[63]. Therefore, this level scheme allows to change the 
sign of Jj_ relative to the {|0), |1)} scheme. The resulting 
values for V, W, J z , and J± are listed in row #3 of Table 
II. 

To calculate the coefficients V, W, J z , and J± in 
the {\m ), Inn)} = {|0), Va|l) + ^/t=a\2)} scheme [Fig. 
3(b); #4 in Table II], we again ignore off-resonant terms. 
Projecting d^d^ onto states |0), |1), and |2), we obtain 

d$p($P « (|0)i<0|ipo + |X>*<1|*a*i + |2)i(2|iM2) 

x(|0) j <0| jM o + |l) J -<l| J -Mi + |2) j <2|,-M2) 

+(/4|01)<10| + /4|02)(20| + A*f a |12><21| + h.c). (37) 

Projecting this on states |0) and |mi), we arrive at 

d )d ] + W+ di - +d i l ) d ( f ) ) = ^B p \m p m p )(m p m p \ 

v 

y A p A q \m p m q ) (m p m q \ + -y-flmomi) (raim | +h.c.$38) 
p,i 

where p,q £ {0, 1} and where the values of A p , B p , and 
J± are listed in row #4 of Table II. Although in the 
present configuration, h(d+d_ +d®d+) does not con- 
tribute and B = 0, we wrote Eq. (38) in this more gen- 
eral form to be able to describe all other configurations 
below. A p can be thought of as an effective dipole mo- 
ment of state \m p ), while B p gives the contribution to the 
interaction from the transition dipole moments between 
the rotor states that make up \m p ). Comparing Eq. (29) 
to Eq. (5), we can read off V = [(A + A^ 2 + B + 
W = (A 2 + B -A 2 1 -B 1 )/2, J z = (Ao-Atf + Bo + Br. 
These expressions hold generally and are listed at the 
top of Table II. The two tuning parameters (dE/B and 
a) can be used, for example, to set W = and J z = J±. 
In particular, at (dE/B, a) = (1.25, 0.74), we get W = 0, 
J z = J± = 0.36d 2 and V = O.U z . Setting W = and 
J z = J± brings Eq. (5) into a form similar to the SU(2)- 
symmetric t-J-V model [58] extended to long-range in- 
teractions. Moreover, as we have noted in Sec. II A 2, the 



value of V = 0.1 Jj is expected to result in a suppression 
of phase separation relative to the original t-J model, in 
which V = -J 2 /4 [8]. 

To find expressions for V, W, J z , and J± in the con- 
figuration {|m >, K)} = {ya|0) + v / T r ^|T>, |1» [Fig. 
3(d); #5 of Table II], we project the resonant terms of 

to obtain 



onto states |1), |0), and |1) 



l 



d^d^) 



(|l) 4 (l|^i + |0) l (0| lM o + | 
x(|l) J (l|^ 1 + |0) 3 (0| jMo 



(II 



,2 

"I IT) 



ill 



(39) 



+ (^oi|01)(10|-^ I |01)(10|- A |- 11 

We note that terms that do not conserve the total M still 
do not contribute since they are all highly off-resonant 
(energy non-conserving). In particular, this assumes, 
that the microwave field is strong enough that |mo) is 

not resonant with i) (or that |<^>i t _i) is shifted away 

using a separate microwave field coupling it, for exam- 
ple, to |02,-2))- Limiting the internal states of the two 
molecules to |mo) and |1), we arrive at Eq. (38) with 
the values of A p , B pi and J± listed in row #5 of Table 
II. The minus signs featured in the expressions for J± 
and J z (when compared to the {|0),-y/a|l) + y/l — a|2)} 
configuration - #4 in Table II) allow to set J± = 
at any dE/B, as well as set J z and/or J± to be neg- 
ative. In particular, if one writes J z — iJlcos^ and 
J± = |J|sin^>, then the ability to achieve any value of 
if) would imply full controllability over J z and Ji and, 
hence, over Eq. (4). And indeed, in a similar configura- 
tion {\m ),\mi)} - {|3),yS|l) + Vl^\2)} [Fig. 3(c); 
#6 in Table II], by tuning a and E, one can achieve 
any value of ip. In particular, in the plane defined by 
dE/B and a, the J z = and J± = lines cross at 
(dE/B, a) = (2.6,0.92), so that all values of ip (and 
hence all four combinations of the signs of J z and J±) can 
be achieved just by going around that point in a circle. 
While this proves that any value of ip can be achieved, 
the resulting values of |J| could be rather small; how- 
ever, it is important to emphasize that for any desired 
if), there is almost certainly a different level configuration 
that gives a larger |J|. 

To achieve an even larger degree of control, one can 
apply two microwave fields. For example [Fig. 3(e); #7 
in Table II], one microwave field can be used to create 
a dressed state |mo) = a\0) + vl — a|l), while another 
microwave field can be used to create a dressed state 
|mi) = 6|1) + vl — b\4>2.-i). Following the same pro- 
cedure as for other level schemes, we arrive at the ex- 
pressions for A p , B p and J± listed in row #7 of Table 
II. In particular, with (dE/B,a,b) = (1.7,0.33,0.81), we 
obtain W = 0, J z = J x = -4V = 0.089d 2 . As we have 
noted in Sec. II A 2, these values of W, V, J z , and J± 
make our model very similar to the original t-J model, 
except the interactions are long-range. Other configu- 
rations with two microwave fields can, of course, also 
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be used to obtain other interesting combinations of co- 
efficients or, possibly, to increase the overall interaction 
strength relative to this example. 

Finally, full controllability can be achieved with three 
microwave fields. In that case, we will have four con- 
trol knobs (three microwave fields and the magnitude of 
the applied DC electric field), which can allow for the 
full control over the four constants V, W, J z , and Jj_. 
In particular, consider the example [Fig. 3(f); #8 in Ta- 
ble II] where the two dressed states are |mo) = \/a\2) + 
V~ a\2) and |mi) = y/b\l) + y/c\T) + y/l - b-c\2). To 
achieve controllability over a, we can apply a o~ field on 
the |2) - |2) transition. At dE/B = 2.97 (see below), 
this transition has frequency 0.3B and a sizable transi- 
tion dipole moment ^2 = — 0.13d. Alternatively, one 
can use a Raman pair of microwaves to couple |2) and 
1 2) via |3) = 1^3,2), in which case the transition dipole 
moments are stronger (/xgg = 0.37<i and /zgg — 0.53(f) and 
the transition frequencies arc larger (~ 6B). To achieve 
controllability over b and c, we can apply a o~ field on the 
|1) — |1) transition (or on the |1) — |2) transition) and a 7r 
field on the |1) — |2) transition. In this configuration, we 
find that the four manifolds V = 0, W — 0, J z — 0, and 
Jl = all intersect at dE/B = (2.97,0.059,0.56,0.38). 
Specifically, to set V = W = J z = J± = 0, it is sufficient 
to set Ai = B\ = Aq + Bq = Jl = 0, which is the pro- 
cedure we followed. Therefore, in a small sphere in the 
4-dimensional (dE/B,a,b,c) space around the intersec- 
tion point of the four manifolds V = 0, W = 0, J z — 0, 
and Jl = 0, one can achieve any value of V, W, J z , 
and J_l up to an overall positive prefactor. While this 
example proves full controllability, the actual magnitude 
of the interaction could be small in this case; however, it 
is important to emphasize that for any desired relation- 
ship between V, W, J Zl and Jj_, there is almost certainly 
a different level configuration and a different choice of 
microwave fields that gives stronger interactions. 

The examples presented here (Table II) are just a very 
small fraction of what is possible. In particular, we would 
like to emphasize that even for the relationships of V, W, 
J z , and Jj_ that we consider, configurations other than 
the ones we present can likely be used to achieve a larger 
overall interaction strength. Similarly, the search for the 
optimal configuration for any given experimental labora- 
tory can take into account the laboratory's constraints 
on the strength of the DC field, on the microwave inten- 
sity, and on the range of available microwave frequencies. 
When designing a configuration to achieve some desired 
relationship between V, W, J z , and J±, various caveats 
can be followed to streamline the search. As one sim- 
ple example of such a caveat, one can ensure that some 
transition dipole moments vanish exactly by using states 
whose N z eigenvalues differ by more than one. This way, 
one can, for example, set J± = independently of the 
strength of the applied DC field. 



VI. STABILITY AGAINST CHEMICAL 
REACTIONS 



We now turn to the discussion of stability of our system 
against chemical reactions. For some species of diatomic 
polar molecules, two absolute (electronic, vibrational, ro- 
tational, hyperfine) ground state molecules cannot react 
to form homonuclear dimers [103, 104]. In that case, one 
may be able to remove the hard-core constraint and con- 
sider Hamiltonians with finite elastic on-site interaction 
(see e.g. Ref. [25, 51]). However, even for these molecules, 
excited states might react [103]. Moreover, the currently 
available molecules, KRb and LiCs, both have exother- 
mic reactions to form homonuclear dimers. Therefore, in 
order to avoid these chemical reactions, it is important 
to ensure that two molecules never sit on the same site. 
There are several ways to prevent two molecules from 
sitting on the same site. First, one can rely on strong 
dipole-dipole repulsion. Specifically, if we use a 2D ge- 
ometry in the X-Y plane or a 3D geometry with Z tun- 
neling shut off, and if we further suppose that the electric 
field direction z is near Z, then, at least for the ground 
rotational state, dipole-dipole repulsion can play the role 
of a hard-core constraint for molecules when they hop in 
the X-Y plane [62] . One may expect that this repulsion- 
induced stability also applies to some situations where 
two rotational states are populated. We will discuss this 
possibility below. Second, sufficiently strong attraction 
between two molecules that sit on the same site should 
also be able to prevent two molecules from hopping onto 
the same site by energy conservation, similar to the ex- 
periments on repulsively bound pairs [105]. Finally, if re- 
action rates [103] are really large, one can also try relying 
on the quantum Zeno effect to provide the hard-core con- 
straint [106, 107]. Therefore, if strong attraction and/or 
the quantum Zeno effect are sufficient to provide stability 
(i.e. strong repulsion is not necessary), our models can be 
extended to the full 3D geometry with tunneling allowed 
along all three directions. 

Let us make an estimate for the suppression of 
chemical reactions caused by the quantum Zeno effect. 
Let w(X) be the ID Wannier function for the poten- 
tial F sin 2 (ifX), where K = 2tt/A and A = 1064 
nm [13]. We can then compute the tunneling ampli- 
tude t=-J dXw(X) [-a^jfsr + V Q sm 2 {KX) w(X - 
A/2) and the on-site chemical reaction rate T — 
[J dXw 4 (X)\ [108], where we take the 3D loss rate 
k 3D = 2x KT 10 cm 3 /s from Fig. 2B of Ref. [23] (which 
is of the same order of magnitude as the theoretical pre- 
dictions of Ref. [109]). We plot t and T in Fig. 6 as 
a function of Vq/Er, where we used the recoil energy 
E R = h 2 K 2 /{2M m ) (2tt)1.4 kHz for mass M m of KRb 
(recall that h = 1). We see that as we increase Vq/Er 
from 5 to 30, r/27r grows from 900 Hz to 5 kHz, while 
t/2-n: drops from 90 Hz to 0.6 Hz. Therefore, t < T is 
satisfied for all the values of Vq considered, and we can 
compute the effective loss rate t 2 /T [108], with the re- 
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suit shown in Fig. 6. We see that even at Vb = 5Er, 
t 2 /r ss (27r)9 Hz, which is already sufficiently slow to 
allow for an experiment to be carried out. The effective 
loss rate falls rapidly to even lower values as we increase 
Vo/Er dropping below 1 mHz at Vq/Er = 30. In partic- 
ular, this means that the simplest {|0), |1)} configuration 
at small dE/B, which gives rise to the J z = V = W = 
model studied in Ref. [7], should be stabilized by the 
quantum Zeno effect, despite the fact that it is not stabi- 
lized by repulsive dipole-dipole interactions (see below). 

It is also important to verify that two molecules on 
neighboring sites would not decay directly due to the 
overlap of their Wannier functions. To do this, we com- 
pute the nearest-neighbor chemical reaction rate = 
k 3D [J dXw 4 (X)] 2 J dXw 2 (X)w 2 {X ~ A/2). As we can 
see from Fig. 6, as we increase Vq/Er from 5 to 30, 
r2/27r drops from 2 Hz to 0.1 mHz, making it negligi- 
bly small. We also see that up to Vq/Er « 30, t 2 /T is 
larger than T2 and, thus, determines the total loss rate. 
By analogy with the interaction-assisted tunneling dis- 
cussed in Appendix B, we can also compute the quan- 
tity T 3 = -k 3D [J dXw 4 (X)] 2 J dXw 3 (X)w(X - A/2), 
which can be thought of as the imaginary part of the tun- 
neling amplitude between two occupied sites. As we can 
see from Fig. 6, T 3 is smaller than t, and, in particular, 
much smaller than T. Therefore, we expect the T 3 pro- 
cess to be suppressed in a way similar to the suppression 
of tunneling t between two occupied sites. 

Let us also make a rough estimate for the strength 
of dipole-dipole interactions for two molecules confined 
to a single site. Taking the dipole moment d of KRb 
and a typical distance of 50 nm between two molecules 
confined to the same site, we get an interaction energy 
E- mt = d 2 /( 47r£ o(50 nm) 3 ) - (2tt)400 kHz, which is much 
larger than the tunneling amplitude t shown in Fig. 6. In 
fact, this interaction energy is even larger than the on- 
site chemical reaction rate T. Therefore, strong dipole- 




FIG. 6. (color online). The tunneling amplitude t, the on-site 
chemical reaction rate F, the effective loss rate t 2 /F, nearest- 
neighbor chemical reaction rate T2, and imaginary part Vs 
of the tunneling amplitude between two occupied sites as a 
function of Vo/Er, where Er is the recoil energy and where 
Vo is the amplitude of the lattice. The vertical axis is in Hz. 
We use A = 1064 nm. 



dipole interactions may further suppress the tunneling 
of molecules between two occupied sites, thus, further 
reducing the loss due to chemical reactions. However, 
a more elaborate calculation [103], which is beyond the 
scope of this paper, is required to fully understand this 
effect. 

Although the quantum Zeno effect and dipole-dipole 
attraction may allow to stabilize the system (as we have 
just described), let us, nevertheless, estimate the stabil- 
ity conditions, assuming we want to rely solely on strong 
repulsive dipole-dipole interactions. Since each molecule 
can be in one of two rotational states, we must ensure 
repulsion for any two-molecule internal state, which will 
significantly restrict the range of parameters at which 
stability is achieved purely by repulsive interactions. Let 
us begin by considering the simplest {|0), |1)} configu- 
ration (#1 in Table II). The terms in Eq. (31) with 
ji = j2 = ja = 3a are the ones that give rise to the 
hardcore constraint. To ensure dipole-dipole repulsion 
between two molecules independently of their internal 
state, two conditions should be satisfied. First, the angle 
Oo that the DC electric field makes with the Z-axis must 
be smaller than sin _1 (l/v / 3) to ensure that Vdd(R) > 
for any vector R in the X-Y plane [see Fig. 1]. Sec- 
ond, we have to require that the term in square brack- 
ets in Eq. (32) is positive for any two-molecule internal 
state. This requirement reduces to a single condition: 
/ioMi > nv wnere we have assumed J± = 2//^ (i.e. 
no tensor shifts). Physically, this condition ensures that 
the two-molecule sing let state (|0)|1) - |l)|0))/-y/2 has 
positive energy. The same analysis can be done for the 
{|1),|3)} configuration (#2 in Table II) and yields the 
stability condition Hi/j, 3 > [i 2 3 . In Fig. 7(a), we show 
the permanent and transition dipole moments that play 
a role in these two configurations; and in Fig. 7(b), we 




FIG. 7. (color online), (a) Permanent and transition dipole 
moments in units of d. (b) Stability curves in units of d 2 . The 
system is stabilized via repulsive dipole-dipole interactions 
when the plotted quantity is positive. 
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plot the stability curves /xqMi — A*oi an d MiA^ - Mi3; whose 
positive values show the regions of stability. We see that, 
in the {|0),|1)} configuration, stability is achieved for 
dE/B > 6, while, in the {|1),|3)} configuration, it is 
achieved for < dE/B < 3.9. The condition dE/B > 6 
requires large electric fields [E > 24(12) kV/cm for KRb 
(LiCs)]. Therefore, it may be easier experimentally to 
achieve stability in the {|1) , |3)} configuration than in the 
{|0), |1)} configuration. The two features of the {|1), |3)} 
configuration that allow it to be stable at small DC elec- 
tric fields are: (1) the fact that /ii and point in the 
same direction at small DC fields and (2) the fact that 
1^13 = for E = 0. We also note that the use of tensor 
shifts to reduce J± may allow one to extend the stability 
range to lower dE/B for some configurations, such as the 
{|0), |1)} configuration. We also point out that the sta- 
bility range for the {|1),|3)} configuration conveniently 
includes the "magic" point for these two states in Fig. 4 
(dE/B 1=3 1.7). Finally, we note that the analysis in the 
present Section can be readily extended to the case when 
microwave fields are applied. 



VII. CONCLUSION 

We derived the t-J-V-W model that governs the be- 
havior of polar alkali dimers in an optical lattice. In 
particular, we showed how microwave fields can be used 
to make the coefficients of the Hamiltonian fully tunable. 
We also described how nuclear spins and the associated 
hyperfine interactions can be used to endow the model 
with another highly controllable (orbital) degree of free- 
dom. The peculiar and highly tunable features of the 
model, such as long-range anisotropic interactions and 
the hyperfine interactions with the nuclear spin, should 
make the system an invaluable resource for studying ex- 
otic manybody phenomena and for providing insights 
into strongly correlated condensed matter systems. 

One feature of the models considered in the present 
manuscript is that, for two nearest-neighbor molecules in 
an optical lattice with 0.5/im spacing [13], dipole-dipole 
interactions are relatively weak (0.4 kHz for KRb and 37 
kHz for LiCs). It would, thus, be convenient to bring 
the molecules closer. First, although the structure of 
molecules is more complicated than that of atoms, and 
inelastic photon scattering rate could vary drastically as 
one tunes the wavelength of the lattice laser [110], we 
believe that lattice spacing down to 200-300 nm will be 
possible. This would increase the dipole-dipole interac- 
tion strength by an order of magnitude. Another promis- 
ing way to achieve closely spaced molecules is to consider 
molecular Wigner crystals [55, 62], which will be the sub- 
ject of future studies. 

Several other extensions of the present work may be 
particularly fruitful. For example, it is straightforward 
to extend the Hamiltonian in Eq. (1) to more than two 
dressed rotational states and, thus, emulate spin S > 
1/2. One can also consider level configurations, in which 



dipole-dipole interaction terms that do not conserve the 
total N z of the two interacting molecules contribute [p 7^ 
in Eq. (A25)] and generate a larger variety of angular- 
dependences in the interaction than that present in Vdd 
[Eq. (2)]. Interaction terms of the form S^S^ can also 
be generated if one uses degenerate dressed states |m ) 
and \mi) allowing one to access, for example, spin models 
beyond the XXZ model. 

Furthermore, while this manuscript is mainly focused 
on quantum simulation applications of the system, ap- 
plications to quantum computation - particularly in the 
context of storing quantum information in the nuclear 
spins [55, 91-93] - can be readily envisioned. By anal- 
ogy with similar proposals for alkaline-earth atoms [88] , 
alkali dimers in an optical lattice may be used, for exam- 
ple, to generate many-body entangled states with appli- 
cations to precision measurements and to measurement- 
based quantum computation. Finally, as another possible 
extension of the present work, we expect that, by anal- 
ogy with Ref . [32] , which treats Rydberg atoms and polar 
molecules on equal footing, our ideas should be extend- 
able to Rydberg atoms. 
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Appendix A: Matrix elements 

In this Appendix, we first present some formulas that 
are useful for evaluating the internal structure of the 
molecules, their interaction with optical and microwave 
fields, as well as their dipole-dipole interaction with each 
other. We then use these formulas to evaluate matrix el- 
ements of the internal molecular Hamiltonian, as well as 
of the dipole-dipole interaction between two molecules. 

Closely following Ref. [94] for most of this Appendix, 
let Tpi(Ai) be a tensor of rank k\ with components p\ 
which operates on angular momentum Ji. Similarly, let 
Tp 2 2 (A2) be a tensor of rank ki with components pi which 
operates on angular momentum J 2 . We assume that Ji 
and J2 commute. We can define the tensor product of 
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T fcl (Ai) and T fe2 (A 2 ) as 

T*(A X , A 2 ) = r* 1 (Ai)Tpl pi (A 2 ) (2fc + I) 1 / 2 



pi 

fci fc 2 fc 
Pi P-Pi ~P 



(Al) 



where the 2x3 matrix in parentheses is the 3-j symbol. 
For k\ = fc 2 = k, we can also define the scalar product of 
T fcl (Ai) and T fe2 (A 2 ) as 

T fc (Ai) • T fc (A 2 ) = ^(-lfT p fe (A 1 )T^(A 2 ). (A2) 



Notice that for k = 1, the spherical and Cartesian scalar 
products agree: T x (Ai) • T X (A 2 ) = A x • A 2 . 

If Ji and J2 couple to form J, we have the follow- 
ing formulas for the reduced matrix elements (reduced 
matrix elements use symbol || instead of | and have no 
dependence on the component indices such as p\, p 2 , and 
p): 



(J 1 ,J 2 ,J\\T^(A 1 )\\J{,J 2 ,J') = 5j 2iJi (-l 



,.7' + .7i + fei + J 2 



Xv /(2J+l)(2J' + l)| J j J Jx £ j(J 1 ||T fc 1 (A 1 )||J0(A3) 

(Ji, J 2 , J\\T k *(A 2 )\\J[, J' 2 , J 1 ) = 5j uJ ,(-l) J+J i +k s +J s 

x v /(2J+l)(2J' + l)|") ■£ £ |(J 2 || T fc2 (A 2 )||J 2 )(A4) 

(Ji,J2, J\\T k (A x , A 2 )\\J[, J 2 , J') = 

f J J' fc 

x v /(2J+l)(2J' + l)(2A: + l)^ J x J[ k x 

{ J2 J' 2 k 2 

x^HT^AOIlJO^HT^Aa)!!^). (A5) 

Here the 2x3 matrix in curly braces is the 6-j symbol, 
and the 3x3 matrix in curly braces is the 9-j symbol. 

The Wigner-Eckart theorem allows to compute matrix 
elements of T*(A), operating on angular momentum J, 
in terms of reduced matrix elements: 



(J,M\T*(A)\J',M') = (-1) J - M ( _ J M k p J M , 
x(J||T fc (A)||J')- 



(A6) 



Three particularly useful sets of reduced matrix ele- 
ments are 

(J\\T\3)\\J') = ^ [J(J + 1)(2 J + 1)] 1 / 2 , (A7) 
J(2J- 1) / J 2 J 



(J||T 2 (J,J)||J'> =5. 



J, J' 



(A8) 



n/6 \-J J 
for any angular momentum J (in particular, for N) and 
(N\\T k {C)\\N') = {-1) N [(2N + l)(2N' + l)] 1 / 2 

J' ( AQ ) 
where T k {C) = C%(6,<t>). 



We will now use Eqs. (A1-A9) to evaluate matrix ele- 
ments of the internal molecular Hamiltonian, as well as 
of the dipole-dipole interaction between two molecules. 
In Ref. [17], which we follow together with Ref. [94] to 
compute the matrix elements, three kinds of bases are 
used. Since we work in the regime where it is sufficient 
to take the expectation value of H^ in a given eigenstate 
of Hq, we will use only two basis sets, as in Ref. [21]: 

\NMM X M 2 ) (uncoupled), (A10) 
\NMIMi) (coupled). (All) 

In both bases, the rotor state \NM) is decoupled from the 
nuclear spin states \I\M-\) and \I 2 M 2 ). The coupled basis 
couples the two nuclear spins and uses \IM/), where I = 
I1+I2, while in the uncoupled basis M\ and M 2 magnetic 
quantum numbers of the two nuclear spins are used. We 
will use the two bases whenever the operator that is being 
considered acts on the nuclear spins. Otherwise - if the 
operator acts only on the rotor degree of freedom - we 
will simply use the basis \NM). 

We begin by computing the matrix elements of Hq [Eq. 
(8)]. Noting that N 2 acts only on the rotor degree of 
freedom, we have 

(NM\N 2 \N'M') = 6 NN ,5 M M'N(N + 1). (A12) 

To evaluate the matrix elements of do, we note that 
dp = Tp(d) = e p ■ d — dCp(0,4>) for all 3 values of 
p = 0, ±1. Thus, for evaluating the DC Stark shift —d E 
and the dipole-dipole interaction between two molecules, 
we need the matrix elements of Cp(9,4>). For evaluat- 
ing the quadrupole hyperfine interaction and the ten- 
sor hyperfine interaction, we need the matrix elements 
of Cp(9,4>). Let us, thus, evaluate the matrix elements 
of Cp(9,4>) for general k. Using Eqs. (A6,A9), we have 

{NM\C>;{e,$)\N'M') = (-1) M [(2N + 1)(2N' + 1)}^ 2 

(A13) 



We now compute the matrix elements of f/hf [Eq. (9)]. 
We begin with Hq. Using the form of Hq in Eq. (9), 
the definitions of T 2 (VE,) and T 2 (Q;), and Eq. (A2), 
we have 

H Q = E(-l) P C 2 (0,0)^|^T 2 p (L i ,I i XA14) 



N' \ 




k 


N' \ 


M' ) 







j 



p. I 



where i sums over the two nuclei. Since we have already 
evaluated the matrix elements of Cp(9,4>), it remans to 
list the matrix elements of T 2 (Ii, Ij). Using Eqs. (A6,A8), 
in the uncoupled basis, they are 



(M i \T 2 v {lul i )m = Ii{21 ^ l) 



x(-l) 



h-M: 



I, 2 Ii \ / Ii 2 1, 
-Mi p Ml ) { -I % U 



(A15) 
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Using Eqs. (A6,A3,A4,A8), in the coupled basis, they are 
(IMjlT^hJ^I'M'j) = 

x (_l)/-Mi+Ji+*.[(2/ + l)(2/' + l)]V2 Wpll 

V6 

7 2 I' \ f L 2 7, 



(A16) 



-Afj p Mj J \ -Ii 7, 

ii i' ij\i(-iy' if *•= i 

if i = 2 ' 



J 4 ^ I ; (-i) r 

\ I I, 2 J t (-1)' 



where j = 2(1) if « = 1(2). 

To compute the matrix elements of 77/ jv , it is sufficient 
[by Eq. (A2)] to compute the matrix elements of N and 
Ii. Using Eqs. (A6,A7), we find that the matrix elements 
of N are 



(NM\T^(N)\N'M') = S NtW (-l) N - M 
x[N(N+ 1)(27V + 1)] 1/2 . 



N 1 A' 
—M p M' 

(A17) 



The matrix elements of Ij in the uncoupled basis are 
[using Eqs. (A6,A7)] 



h 1 I'i 

Mi p M[ 



x[7 l (7 l + l)(27, + l)] 1 / 2 . (A18) 
In the coupled basis, they are [using Eqs. (A6,A3,A4,A7)] 

7 1 V 



\/-M/+/i+/2 



(7M / |T;(I i )|7'A7|) = -(-l) p M , 

x[(27 + l)(27' + l)] 1 /2|^ £ h J[J 4 (/ 4 + l)(2J 4 + l)]i 



if z = 1 



(A19) 



xi ^ 

X l (-1)' if i = 2 

where j = 2(1) if i = 1(2). 

We now turn to 77 t . Using Eq. (A2), 

77 t = -c 3 y6^(-l)fC 2 p (0,0)T p 2 (I 1 ,I 2 ). (A20) 



Thus, since we have already evaluated the matrix ele- 
ments of Cp, it remains to evaluate the matrix elements 
of Tp (Iijla). In the uncoupled basis, they are [using Eqs. 
(A1.A18)] 

( Mi M 2 1 T p 2 (Ii,I 2 ) I M[M' 2 ) = (-\)h-M 1+ i 2 -M 2 - v 



x 



[57 1 (7 1 + l)(27j + 1)7 2 (7 2 + 1)(27 2 + I)] 1 ' 2 



x E 

Pi=— l 

' h 1 7 2 



Pi p-pi -p J V -Mi pi A7( 



x 



-M 2 p-p x M 2 



(A21) 



In the coupled basis, they are [using Eqs. (A6,A5,A7)] 

(7M / |r 2 (I 1 ,I 2 )|7'M;) = (-1) J " M ' [5(27 + 1) (27' + l)]i 
x[/i(7i + l)(27 a + 1)7 2 (7 2 + 1)(27 2 + I)] 1 ' 2 



(A22) 



Finally, the matrix elements of 77 sc in the uncoupled 
basis are [using Eqs. (A2,A18)] 

(Mi M 2 1 77 sc | M[M' 2 ) = Ci (-\)h-M 1+ i 2 -M 2 
x[h{h + l)(27i + 1)7 2 (7 2 + 1)(27 2 + 1)]5 

7 2 1 I 



E 



0=-l 



h 1 h 

-Mi p M( 



-M, 



p M> J (A23) 



In the coupled basis, they are (|7i — 7 2 | < 7 < 7i + 7 2 is 
assumed) 



(IMj\H sc \I'Mi) = c 4 5 ir S MlM , 

x ^[7(7 + 1) - 7i(7i + 1) - 7 2 (7 2 + 1)]. 



(A24) 



We conclude this Appendix by presenting a conve- 
nient expression for dipole-dipole interaction between 
molecules 1 and 2 separated by R = (R,9',<f)'), where 
9' and </>' are the spherical angles of R in the x-y-z coor- 
dinate system, which is defined with respect to the direc- 
tion z of the applied DC electric field. This expression 
is: 



77, 



V6 



dd 



47re 7? 3 
~4ire R 3 



T 2 (C)-T 2 (dW d< 2 )) 

2 

]T (-l)^ p (C)T p 2 (d«,d(%25) 

p=-2 



^Yk^O'A') and 



where T 2 (C) = C* (#,<!/) = 

where we used Eq. (A2). In the present manuscript, we 
only use the p = component, for which [using Eq. (Al)] 

T 2 (d« d (2) ) = j= (d^df + 2^ 1) 4 2) + d^d^6) 



T 2 (C) = -(3cos 2 0'-l). 



(A27) 



It is easy to see from Fig. 1 that cos 9' = R • z 
sin O cos($ - 4> ). 



Appendix B: Interaction- Assisted Tunneling 

In this Appendix, we analyze the small corrections to 
the two approximations made to arrive at Eq. (32). The 
two approximations were: (1) the extent of the Wannier 
function w is much smaller than the distance between 
the sites, and (2) only terms involving two sites and con- 
serving the number of molecules on each site contribute. 
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Making the second approximation, corrections to the first 
approximation lead to the replacement of Vdd(Ri — Rj) 
in Eq. (32), when i and j are close to each other, with a 
more complicated dependence on i, j and on the internal 
state of the molecules at i and j. However, even the sep- 
aration of nearest neighbor sites (~ 500 nm) is typically 
much larger than the extent of Wannier functions in a 
deep lattice (~ 50 nm). 

Corrections to the second approximation result in 
interaction-assisted tunneling. Interaction-assisted tun- 
neling corresponds to the terms in Eq. (31) where ji = j\, 
while j2 and j 3 are nearest neighbors not equal to j\. 
Terms where j'2 and j 3 are not nearest neighbors are 
much smaller. Terms where ji = j 3 while ji and j'4 
are nearest neighbors are identical to the case we are de- 
scribing, giving an overall factor of 2. Another factor of 
2 comes from the fact that (j'2,.73) counts each nearest- 
neighbor pair only once. The resulting Hamiltonian is 
approximately equal to (the subscript in iJj at stands for 
interaction-assisted tunneling) 

Hmt = 2 ^ lull VxUuh, h ) S f n c] 2la , c hQa , + h. 

ji<r' \ 

(h, J3>#ji 

h (Bl) 

mm' ) 



where 



J d 3 Rd 3 R'Vdd(R - R') 
xwjio(R)w Ji i(R)w J2 i(R')w J3 o(R'), 
J d 3 Rd 3 R'Vd<j(R - R') 



X!0 



! m (RK ffl( (RV(R'). (B2) 



Physically, the interaction-assisted tunneling means that 
the presence of a molecule on site j\ assists in a "tun- 
neling" of a molecule from site j 3 to site j'2 ■ In the term 
proportional to V mm t, this "tunneling" does not change 
the internal states of the two molecules, while in the term 
proportional to Vj_, it is accompanied by an exchange of 
a rotational excitation between the two molecules. 

Let us compare the magnitude of the tunneling ampli- 
tude t to the magnitude of the interaction-assisted tun- 
neling. In the deep-lattice limit, the tunneling amplitude 
t in a ID potential Vq sin 2 (if X) (where K = 2ir/X), is re- 
duced relative to the recoil energy Er — h 2 K 2 /(2M m ) by 

a factor porportional to exp (^—2^/Vq/ Ejij [111]. Since 



interaction-assisted tunneling also involves an overlap of 
Wannier functions on neighboring sites, we may expect it 
to fall off similarly with increasing Vq. At the same time, 
the reference energy scale for interaction-assisted tunnel- 
ing is i?dd = ^ 2 /(47reo(A/2) 3 ), the strength of dipolc- 
dipole interaction between nearest-neighbor sites. For 
KRb with A = 1064 nm [13], E R w (2tt)1.4 kHz > E dd w 
(2-7r)0.3 kHz, so we may expect the interaction-assisted 
tunneling to be smaller than the usual tunneling. 

To be more precise, in Fig. 8, we compare the magni- 
tude of interaction-assisted tunneling to the usual tun- 
neling. Let w(X) be the ID Wannier function for 
the potential Vq sin 2 (KX). For the case when the 
(X, Y) coordinates of the three sites are (in units of 
a = A/2) j, = (0,0), h = (1,0), and j 3 = (2,0), 
we estimate the amplitude of interaction-assisted tun- 



dX^,w{X - a)w(X - 2a). 



neling as V x = - f a °j!^ 
We do not integrate from X = —00 to avoid inte- 
grating over the singularity of the 1/A 3 potential at 
X = 0, which is unphysical and stems from the fact 
that 1/X 3 interaction breaks down at small X. For 
the case ji = (0,1), j 2 = (0,0), and j 3 = (1,0), we 
estimate the amplitude of interaction assisted tunneling 
as F 2 = JZo dX 4^(^+^)3/2 w(X)w(X ~ a). Solving 
numerically for w(X) and for the tunneling amplitude 
t, assuming the mass and dipole moment of KRb and 
A = 1064 nm, in Fig. 8, we plot V\, V 2 , and t as a 
function of Vq/Er. We see that for the presented val- 
ues of Vo/Er, the interaction- assisted tunneling is at 
least 10 times weaker than the usual tunneling ampli- 
tude t, which confirms our expectations and allows to 
ignore interaction-assisted tunneling. 




V /E. 



R 



FIG. 8. (color online). The tunneling amplitude t and 
interaction-assisted tunneling amplitudes Vi and V2 as a func- 
tion of Vo/ En, where En is the recoil energy and where Vb is 
the amplitude of the lattice. The vertical axis is in Hz. We 
use A = 1064 nm and the mass and permanent dipole moment 
of KRb. 
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